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Chapter  1 

Executive  Summary 


The  project  entitled  “/I  Computational-experimental  Approach  to  Hierarchical  Modeling  of  Damage 
and  Failure  in  Non-uniform  Composite  Materials”  began  in  February  1998.  A  no-cost  extension 
was  granted  to  continue  this  project  till  June  30,  2001.  During  the  period  of  this  grant,  substantial 
progress  has  been  made  in  advancing  the  state  of  the  art  in  multiple  scale  modeling  and  damage 
by  interfacial  debonding  and  fiber  cracking  in  composite  materials.  Research  has  been  conducted 
in  a  few  distinct  areas  that  are  delineated  below.  A  list  of  publications,  acknowledging  this  grant 
is  provided  in  chapter  2. 

I. Interfacial  Debonding  Analysis  in  Multiple  Fiber  Reinforced  Composites 
(Details  provided  in  Chapter  3) 

Decohesion  at  multiple  fiber  interfaces  of  elastic  fiber  reinforced  composites  is  modeled  by  the 
Voronoi  cell  finite  element  model  (VCFEM)  in  this  chapter.  Interfacial  debonding  is  accommo¬ 
dated  by  cohesive  zone  models,  in  which  normal  and  tangential  springs  tractions  are  expressed  in 
terms  of  interfacial  separation.  Model  simulations  are  compared  with  results  from  experiments, 
performed  using  cruciform  specimens  of  single  and  multiple  fiber  polymer-matrix  composite?.  An 
inverse  problem  is  solved  to  calibrate  the  cohesive  zone  parameters.  Debonding  at  fiber  matrix 
interfaces  aire  simulated  for  different  architectures,  volume  fractions  and  boundary  conditions,  to 
understand  the  influence  of  microstructural  morphology  and  boundary  conditions  on  the  decohe¬ 
sion  process. 

II.  A  Multi-level  Computational  Model  for  Multi-scale  Damage  Analysis  in  Composite 
and  Porous  Materials 

(Details  provided  in  Chapter  4) 

.A.n  adaptive  multi-level  methodology  is  developed  in  this  chapter  to  create  a  hierarchy  of  com¬ 
putational  sub-domains  with  varying  resolution  for  multiple  scale  problems.  It  is  intended  to 
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concurrently  predict  evolution  of  variables  at  the  structural  and  microstructural  scales,  as  well 
as  to  track  the  incidence  and  propagation  of  microstructural  damage  in  composite  and  porous 
materials.  The  microstructural  analysis  is  conducted  with  the  Voronoi  cell  finite  element  model 
(VCFEM),  while  a  conventional  displacement  based  FEM  code  executes  the  macroscopic  analysis. 
The  model  introduces  three  levels  in  the  computational  domain  which  include  macro,  macro-micro 
and  microscopic  analysis.  It  differentiates  between  non-critical  and  critical  regions  and  ranges  from 
macroscopic  computations  using  continuum  constitutive  relations  to  zooming  in  at  hotspots  for 
pure  microscopic  simulations.  Coupling  between  the  scales  in  regions  of  periodic  microstructure 
is  accomplished  through  asymptotic  homogenization.  An  adaptive  process  significantly  increases 
the  efficiency  while  retaining  appropriate  level  of  accuracy  for  each  region.  Numerical  examples 
are  conducted  for  composite  and  porous  materials  with  a  variety  of  microscopic  architectures  to 
demonstrate  the  potential  of  the  model. 

III.  Experimental-Computational  Investigation  Of  Damage  Evolution  In  Discontinu- 
ously  Reinforced  Aluminum  Matrix  Composite 

(Details  provided  in  Chapter  5) 

This  chapter  deals  with  a  combined  experimental-computational  approach  to  study  the  evolution 
of  microscopic  damage  to  cause  failure  in  commercial  SiC  particle  reinforced  DRA’s.  Determi¬ 
nation  of  aspects  of  microstructural  geometry  that  are  most  critical  for  damage  nucleation  and 
evolution  forms  a  motivation  for  this  work.  An  interrupted  testing  technique  is  invoked  where 
the  load  is  halted  in  the  material  instability  zone,  following  necking  but  prior  to  fracture.  Sample 
microstructures  in  the  severely  necked  region  cire  microscopically  examined  in  3D  using  a  serial 
sectioning  method.  The  micrographs  are  then  stacked  sequentially  on  a  computer  to  reconstruct 
3D  microstructures.  Computer  simulated  equivalent  microstructures  with  elliptical  or  ellipsoidal 
particles  and  cracks  are  constructed  for  enhanced  efficiency,  which  are  followed  by  tessellation  into 
meshes  of  2-D  and  3-D  Voronoi  cells.  Various  characterization  functions  of  geometric  parameters 
are  generated  and  sensitivity  analysis  is  conducted  to  explore  the  influence  of  morphological  pa¬ 
rameters  on  damage.  Micromechanical  modeling  of  2D  micrographs  are  conducted  with  Voronoi 
Cell  Finite  Element  Method  (VCFEM).  Inferences  on  the  initiation  and  propagation  of  damage  are 
made  from  the  2D  simulations.  Finally,  the  effect  of  size  and  characteristic  lengths  of  representative 
material  element  (RME)  on  the  extent  of  damage  in  the  model  systems  is  investigated. 
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Chapter  3 


Interfacial  Debonding  Analysis  in 
Multiple  Fiber  Reinforced  Composites 

3.1  Introduction 

Affordability  issues  are  gradually  shifting  design  practice  away  from  the  assembly  of  small  sub¬ 
components  towards  large  unitized  structures.  In  this  paradigm  shift,  greater  emphasis  is  placed 
on  design  based  on  analysis,  rather  than  on  laboratory  tests.  The  micro-mechanisms  of  failure  in 
composite  materials,  e.g.  polymer  or  organic  matrix  composites  are  generally  known.  They  involve 
either  fiber  splitting  or  fiber-matrix  interface  decohesion,  followed  by  matrix  cracking.  The  damage 
mechanisms  are  sensitive  to  local  morphological  parameters  like  volume  fraction,  size,  shape  and 
spatial  distribution  of  reinforcements,  interfacial  strength  and  process-related  defects.  Fibers  and 
interfaces,  in  regions  of  clustering  or  preferential  alignment,  are  subjected  to  increased  local  stresses 
and  have  a  greater  propensity  to  undergo  damage  nucleation  than  those  in  dilute  regions.  A  number 
of  failure  prediction  models  for  composite  laminates  are  phenomenological  and  have  been  based  on 
empirical  methods  or  ply  level  fracture  mechanics.  These  macroscopic  models  are  not  capable  of 
relating  the  failure  process  to  local  interactions  and  stress  concentrations.  Because  of  limited  large- 
scale  tests  in  the  emerging  design  environment,  it  is  critical  that  a  robust  failure  analysis  model 
evolve  from  the  details  of  the  microstructure  and  incorporate  the  micromechanics  of  damage  modes. 

Analysis  of  damage  evolution  by  decohesion  of  fiber-matrix  interfaces  in  multiple  fiber-reinforced 
microstructures  is  the  objective  of  the  present  study.  The  effect  of  weak  bonding  or  debonded  in¬ 
terface  on  the  mechanical  properties  have  been  studied  by  several  investigators  e.g.  [4,  13,  26], 
using  simplified  models  for  representing  imperfect  conditions  through  traction  discontinuities.  The 
propagation  of  interfacial  cracking  or  decohesion  at  fiber- matrix  interfaces  has  been  successfully 
modeled  by  a  number  of  researchers  using  the  cohesive  volumetric  finite  element  methods.  Among 
the  important  contributions  to  the  field  of  damage  evolution  by  normal  and  tangential  separation 
are  those  by  Needleman  [22,  23,  24],  Tvergaard  [29,  30].  Allen  et.  al.  [2,  18],  Lissenden  et.  al. 
[17],  Geubelle  et.  al.  [9,  15],  Walter  et.  al.  [32],  among  others.  A  majority  of  these  studies  have 
used  unit  cell  models,  which  assume  that  the  material  is  constituted  of  periodic  repetition  of  single 
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cells.  Displacement  based  finite  element  analyses,  with  the  inclusion-matrix  interface  represented 
through  traction-displacement  constitutive  models,  are  used  to  predict  the  onset  and  growth  of 
debonding  .  While  these  models  provide  valuable  insights  into  the  microstructural  damage  pro¬ 
cesses,  the  simple  microstructure  are  ineffective  in  addressing  the  interaction  between  fibers,  effects 
of  clustering,  alignment,  relative  sizes  etc.,  that  are  often  observed  in  micrographs.  To  overcome 
this  limitation,  Zhong  and  Knauss  [33]  have  proposed  a  hybrid  discrete-continuum  approach  in 
which  discrete  particle  interactions  with  damage  evolution  axe  modeled,  accounting  for  particle  size 
and  spacing. 

In  this  work,  decohesion  at  multiple  fiber  interfaces  is  modeled  for  polymer  matrix  compos¬ 
ites  by  the  Voronoi  cell  finite  element  model  (VCFEM),  developed  by  Ghosh  and  coworkers  [34]. 
VCFEM  has  been  established  as  an  effective  tool  for  efficient  and  accurate  modeling  of  non-uniform 
microstructures  with  heterogeneities  of  arbitrary  shapes,  sizes  or  dispersions  with  perfect  interfaces 
in  [35,  34,  21].  In  VCFEM,  the  computational  model  evolves  by  tessellation  of  the  microstruc¬ 
ture  to  generate  a  mesh  of  multi-sided  Voronoi  polygons  as  shown  in  figure  3.1.  The  Voronoi  cell 
formulation  is  augmented  in  this  work  to  incorporate  interfacial  debonding  through  the  introduc¬ 
tion  of  cohesive  zone  models.  This  is  accomplished  by  permitting  the  elastic  matrix  and  inclusion 
phases  to  be  connected  at  nodes  by  normal  and  tangential  cohesive  springs.  The  debonding  pro¬ 
cess  is  assumed  to  be  quasi-static,  neglecting  inertia.  The  analysis  is  assumed  to  be  in-plane  and 
post-debonding  friction  is  neglected  in  this  study.  Simulations  with  the  model  aie  compared  with 
results  from  experiments  performed  using  single  and  multiple  fiber  cruciform  specimens.  An  inverse 
problem  is  solved,  using  minimization  techniques  to  calibrate  cohesive  zone  parameters  for  compar¬ 
ison  studies.  Finally,  debonding  at  multi-fiber  interfaces  are  simulated  for  different  architectures, 
volume  fractions  and  boundary  conditions  to  understand  their  effect  on  decohesion. 


3.2  Interfacial  Debonding  with  the  Voronoi  Cell  FEM 

The  Voronoi  cell  finite  element  model  (VCFEM)  for  a  heterogeneous  domain  with  a  dispersion 
of  inclusions  or  voids,  implements  a  mesh  of  Voronoi  polytopes.  The  Voronoi  cells,  surrounding 
each  heterogeneity,  axe  generated  by  a  surface  based  tessellation  algorithm  [11,  14]  to  yield- mul¬ 
tiple  edges  depending  on  the  number  of  neighboring  heterogeneities.  Each  element  in  VCFEM 
consists  of  the  heterogeneity  and  its  neighborhood  region  of  matrix.  No  further  refinement  of  the 
element  is  required  with  this  formulation.  VCFEM  has  been  successfully  developed  by  Ghosh 
et.  al.  [35,  34,  21]  for  accurate  stress  and  deformation  analysis  of  elastic  and  plastic  materials 
with  perfect  interfaces.  It  uses  an  assumed  stress  hybrid  finite  element  formulation  [35,  21]  with 
specially  developed  equilibriated  stress  fields  consistent  with  micromechanics.  VCFEM  is  able  to 
significantly  enhance  computing  efficiency  for  complex  microstructures  compared  to  conventional 
displacement  based  FEM  packages.  Large  regions  of  the  microstructure  axe  thereby  easily  modeled 
by  this  method.  A  brief  account  of  the  extension  of  VCFEM  to  accommodate  theinclusion-matrix 
interfacial  debonding  is  presented  next. 

Consider  a  typical  multiphase  domain  or  representative  volume  element  Ct  consisting  of  N  in- 
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elusions  (see  figure  3.1a),  each  contained  in  a  Voronoi  cell  element  fig,  as  shown  in  figure  3.1b. 
The  matrix  and  the  inclusion  or  reinforcement  phases  in  each  element  cure  denoted  by  Qm  and  Qc 
respectively,  i.e.  fig  —  LJ  ^c*  Each  element  boundary  ^fig  is  assumed  to  be  comprised  of  the 
prescribed  traction  and  displacement  boundaries  Ttni  and  r„m  respectively,  and  the  inter-element 
boundary  Tm,  i-e-  ^fig  =  remljrumUrm.  For  allowing  decohesion  of  the  matrix-inclusion  inter¬ 
face.  an  incompatible  displacement  field  is  facilitated  across  the  interface  through  a  set  of  connected 
node-pairs.  The  nodes  in  each  pair  belong  to  the  matrix  and  inclusion  boundaries  and  dfl^, 
respectively.  has  an  outward  normal  {=«'"),  while  n®  is  the  outward  normal  to  dfle-  It 
should  be  noted  that  the  interfacial  zone  has  zero  thickness  prior  to  deformation,  but  the  nodes 
may  separate  with  progression  of  deformation  and  onset  of  decohesion. 


An  incremental  VCFEM  formulation  accommodates  evolving  interfacial  damage  with  chang¬ 
ing  applied  loads,  deformation  and  stress  fields.  Let  cr"*  and  cr^  be  the  equilibriated  stress  fields 
and  e"*  and  the  corresponding  strain  fields  in  the  matrix  and  inclusion  phases  of  each  Voronoi 
element  respectively.  The  prefix  A  corresponds  to  increments.  Furthermore,  let  u,  and 
denote  kinematically  admissible  displacement  fields  on  dQe-  and  dQ‘  respectively.  A  comple¬ 
mentary  energy  functional  may  be  expressed  for  each  element  in  terms  of  increments  of  stress  and 
boundary/interface  displacement  fields  as: 


ng((T,  A<7,  u,  Au)  =  -  /  A5(^'",  Ao-'”)  dQ-  f  A<t^)  dQ  -  [  e'"  :  Aa 

^  Jam 


dO. 


-L 


:  A(T^  dQ  + 


f  (tr^  +  Ao-'")  •  n"  •  (u"*  +  Au"')  ddQ  -  [  (t  -t-  At)  •  (u”'  +  Au"*) 

JdQe 


dT 


-  f  (o-’"  -h  Act'")  •  n*"  •  (u"'  +  Au'")  ddQ  +  [  +  Acr'^)  ■  n‘=  •  +  Au*") 

-  /  T-dCu-  -  <)  ddQ 

r  r(u7‘  +  Au7'-uf-Auf) 

-  /  T^diuT  -  ut)  ddQ 


ddQ 


(3.1) 


where  B  is  the  complementary  energy  density  and  the  superscripts  m  and  c  correspond  to  variables 
associated  with  the  matrix  and  inclusion  phases.  The  different  terms  in  the  right  hand  side  of 
equation  (4.21)  are  included  to  provide  weak  forms  of  essential  governing  equations  of  the  problem. 
The  fifth  term  corresponds  to  inter-element  traction  reciprocity  while  the  sixth  term  accounts  for 
the  boundary  traction.  The  last  two  terms  provide  the  work  done  by  the  interfacial  tractions 
T”*  =  due  to  interfacial  sepeuration  (u"*  -  u'^).  where  Tn  and  Tj  are  the  normal  and 

tangential  components.  The  traction  on  the  matrix  interface  dQ^  is  related  to  the  stresses  and 
interface  normals  as 

T”*  =  cr"*  •  n"*  =  -o-'"  •  n"  =  -cr"  •  n"  (3.2) 


The  interf3u:ial  response  is  described  by  constitutive  relations  that  prescribe  the  dependence  of 
normal  and  tangential  tractions  on  components  of  interfacial  separation.  The  total  energy  for  the 
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(3.3) 


entire  heterogeneous  domain  is  obtained  by  adding  the  energy  functionals  for  N  elements 

s 

n  =  y:n, 

e=l 

Equating  the  variation  of  lie  in  equation  (4.21)  with  respect  to  stress  increments  A<t"‘  and  Ao-^' 
to  zero,  yields  the  element  displacement  compatibility  relations  in  each  of  the  phases  flm  and  Qc- 
Furthermore,  setting  the  variation  of  FI  in  equation  (3.3)  with  respect  to  the  independent  boundary 
displacements  Au,  Au”*  and  Au'^  to  zero,  yield  the  traction  reciprocity  conditions  on  the  interele¬ 
ment  boundaries  (r„i)  and  traction  boundaries  (Ft^),  and  Hber-matrix  interfaces,  and 
respectively. 


Independent  assumptions  on  stress  increments  Aff  are  made  in  the  matrix  and  reinforcement 
phases  to  accommodate  stress  jumps  across  the  interface.  In  two-dimensional  analysis,  the  Airy  s 
stress  function  #(i,y)  is  used  as  a  convenient  tool  for  deriving  equilibriated  stress  increments.  An 
essential  micromechanics  observation,  that  interfacial  stress  concentrations  depend  on  the  shape  of 
the  inclusion,  has  been  incorporated  in  the  choice  of  stress  functions,  by  Moorthy  and  Ghosh  [34] 
through  the  decomposition  of  the  matrix  stress  functions  into  (a)  a  purely  polynomial  function 
and  (b)  a  reciprocal  function  ■  The  inclusion  stress  functions  are  admitted 

as  polynomial  function  (#^  =  The  pure  polynomial  function  ^pdy  accommodates  the 

far  field  stress  in  the  matrix  and  inclusion  and  are  written  as: 


=  and  ^po(y  =  E 


(3.4) 


P^Q 


The  reciprocal  stress  function  in  the  matrix  facilitates  stress  concentration  at  the  interface 
accounting  for  inclusion  shape  and  decays  at  large  distances  from  the  interface. 


n  1  /  AR^  Aff^  \ 

The  function  f{x,y)-  is  a  specially  mapped  radial  coordinate  that  has  the  properties 


(3.5) 


/(x,y)  =  1  on  and  ^  -^0  as  (x,y)  -»•  oo 


(3.6) 


For  2D  elliptical  heterogeneities,  /(x,  y)  can  be  constructed  by  conformal  mapping  while  for  3D, 
ellipsoidal  coordinates  can  be  chosen  to  represent  this  function.  Stress  increments  in  the  matrbc 
and  inclusion  phases  of  the  Voronoi  cell  elements  may  he  obtained  by  differentiating  the  stress 
functions  with  respsect  to  x  and  y,  to  yield  expressions  of  the  form: 


Ao-  ] 

f  AaJ,  1 

Ao- 

=  [P"*]{A^'”}  and  Ao^^^ 

,  Aa- 

1  J 

[Pl{Ar} 


(3.7) 
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where  [P"*]  and  [P*^]  are  the  stress  interpoltation  matrices  in  the  matrix  and  inclusion.  The  bound¬ 
ary  displacements  are  generated  by  interpolation  in  terms  of  nodal  displacements  on  Sfig,  and 
dCl^  using  conventional  lineeir  or  quadratic  shape  functions. 

{Au}  =  [L^]{Aq}  cmdQe,  {Au"*}  =  [L=]{Aq"*}  on  dQ^  {Au^}  =  [L'^]{Aq=}  on  an^3.8) 

where  {Aq}.  {Aq'"}  and  {Aq*"}  are  the  generalized  nodal  displacement  vectors.  The  stress  and 
displacement  interpolations  are  substituted  in  equations  4.21  and  3.3.  Subsequently  stationarity 
conditions  of  these  equations  axe  evaluated  with  respect  to  the  stress  parameters  A/3  and  A^  ,  ^■tid 
displacement  parameters  {Aq},  {Aq^}  and  {Aq<=}  to  yield  the  stress  and  displacement  solutions 
in  each  element.  The  details  of  VCFEM  formulation  and  solution  methodology  are  presented  in 
[35.  34]  and  is  not  repeated  here. 

3.2.1  Cohesive  Zone  Models  for  Interfacial  Decohesion 

Modeling  the  decohesion  of  matrix-reinforcement  interfaces  is  conveniently  accomplished  using 
crack  initiation  and  propagation  criteria  based  on  the  evaluation  of  traditional  fracture  mechanics 
parameters  K  and  J.  An  alternative  approach  has  been  pioneered  by  Dugdale  [8],  Barenblatt  [3], 
Rice  [27]  and  others  to  avoid  crack  initiation  and  growth  conditions.  This  approach  uses  cohe¬ 
sive  zone  interfacial  models  to  depict  fracture  as  a  phenomenon  of  progressive  separation  across  an 
extended  crack  tip  or  cohesive  zone,  that  is  resisted  by  cohesive  tractions.  Traction  across  the  inter¬ 
face  reaches  a  maximum,  subsequently  decreases  and  eventually  vanishes  with  increasing  mterfacid 
separation,  signaling  complete  decohesion.  Dimensional  considerations  introduce  a  characteristic 
length  in  these  models.  The  interface  mechanical  response  is  specified  in  terms  of  critical  interfacial 
strength  and  work  of  sepauration  per  unit  area  and  no  additional  failure  criteria  are  required.  A 
number  of  cohesive  zone  models  have  been  developed  to  characterize  interfacial  decohesion  at  the 
continuum  scale.  Needleman  [22,  23,  24]  has  proposed  a  potential  based  framework  to  describe 
debonding  initiation  through  complete  separation  in  the  cohesive  zone.  A  similar  potential  based 
model  has  also  been  developed  by  Tvergaard  [29,  30].  Cohesive  zone  models  utilizing  nonlinear 
spring  models  to  depict  the  interfacial  failure  process  have  been  proposed  by  Knauss  and  coworkers 
[31.  33).  Geubelle  et.  al.  [9,  15].  From  thermodynamic  considerations,  Costanzo  and  Allen  [6,  7] 
have  postulated  rate  dependent  cohesive  zone  models  with  internal  state  variables  to  represent  mi¬ 
croscopic  dissipation  mechanisms.  Ortiz  et  al.  [5,  25]  have  developed  a  class  of  irreversible  cohesive 
laws  from  potential  or  free  energy  functions  for  tracking  dynamically  growing  cracks. 

In  this  work  the  rate- independent  cohesive  zone  models  of  Needleman  [22,  23,  24]  and  Geubelle 
et.  al.  [9.  15]  are  incorporated  in  the  Voronoi  cell  FEM  to  model  initiation  and  progressive  debond¬ 
ing  of  matrix-inclusion  interfaces  in  composites.  The  loading  portion  of  the  interfacial  decohesion 
is  manifested  by  these  laws.  In  Needleman’s  model,  the  normal  and  tangential  components  of  inter¬ 
facial  traction  are  derived  from  a  potential,  that  is  expressed  in  terms  of  polynomial  or  exponential 
functions  of  displacement  jumps  or  separation  across  the  interface.  The  potential  using  polynomial 
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(cubic)  functions  is  expressed  in  the  form; 


27 

(j){Un,Ut)  = 


2 


where  Un  and  ut  are  the  normal  and  tangential  components  of  displacement  jump  at  the  interface. 
The  cohesive  zone  parameters  are  cTman  the  maximum  traction  cairried  by  the  interface  undergoing 
a  purely  normal  separation  i.e.  ut  =  0,  (5*,  a  characteristic  length  and  a,  the  ratio  of  interface  shear 
to  normal  stiffness  of  the  interface.  The  tractions  are  obtained  by  dlrterentiating  the  potential  as 


Tn  = 
Tt  = 


d<t)  _  27  f  [i  _  2 4- a 


d(p 

dUn 

d4> 

dui 


-  1 


}  (3.10) 


-  -YCTmaxla  )  1  + 


O' 


For  Un  >  <5*.  the  traction  components  r„  and  Tt  are  zero  in  this  model.  In  an  alternative  form,  the 
potential  is  expressed  using  exponential  functions  of  the  separation  as: 


(/i(Un,U()  —  “CTniax*^  ^  {1 


1  f  0  ^  2  ( 


where  z  =  1|^.  Thus  the  interfacial  traction  components  are 


e[z  ( -  -az^  ( 

dun  ~  ^  \S-J  2^^  U‘V  ^ 

d<p 


(3.11) 


(3.12) 


Tn  =  -■ 

JUn 

Coupling  between  normal  and  tangential  components  is  achieved  through  the  functions  of  Un  and  Ut- 


The  cohesive  failure  model  by  Geubelle  et.  al.  [9.  15]  incorporates  a  bilinear  form  of  the 
traction-displacement  jump  relation  for  mode  mixity  of  applied  tractions.  The  traction-separation 
law  is  expressed  as 


Tn 


Tt 


for(5  <  Smax 

Omax  - 

‘'T  iovS>Smax 


St 


for  S  ^  ^mcLX 
for  (5  >  ^max 


(3.13) 


where  Un  and  ut  are  the  normal  and  tangential  components  of  interfacial  displacement  jump  and 
(5n,  St  and  6  are  the  corresponding  non-dimensional  variables,  defined  as 
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(3.14) 


j.  Un 

C  ' 

yt 

•Sn  =  -7 
K 

The  parameters  and  are  critical  values  of  normal  and  tangential  separation,  for  which  complete 
debonding  is  assumed  to  occur.  The  parameter  cJmai  corresponds  to  the  maximum  value  of  normal 
traction  Tn  for  a  normal  displacement  jump  of  u„  =  SmaxUn-  The  maximum  absolute  value  of  \Tt\ 
is  Tmaz  =  (^max'§  which  occurs  at  |ut|  =  ^^maxul  ■  Since  the  work  of  separation  per  unit  area  on 
the  interface  in  the  normal  and  tangential  directions,  are  expressed  as 


Tn 


2  ^  rnax  ^ 


and 


^  C 


(3.15) 


The  maximum  tangential  traction  expression  implies  that  r„  =  Fj,  or  that  the  critical  energy  re¬ 
lease  rates  are  the  same  for  mode  I  and  mode  II.  This  constraint,  in  a  sense,  creates  an  equivalence 
between  the  potential  based  models  and  the  non-linear  spring  relations.  Parameters  in  both  classes 
of  models  may  be  adjusted,  so  that  they  possess  the  same  essential  features,  e.g.  the  fracture 
energy  or  the  area  under  the  cohesive  law,  the  peak  cohesive  traction  etc..  The  parameters 
or  6’  typically  introduce  a  characteristic  length  of  the  material  microstructure.  A  major  difference 
between  the  two  models  is  that  the  Geubelle’s  model  admits  a  jump  in  the  slope  of  the  traction- 
displacement  curve,  while  the  Needleman's  polynomial  and  exponential  models  have  continuous 
slopes. 


Both  of  these  laws  are  reversible  i.e.,  they  retrace  the  traction-displacement  curve  upon  unload¬ 
ing.  It  has  been  argued  (see  Ortiz  et.  al.  [25,  5])  that  the  decohesion  process  is  expected  to  exhibit 
some  irreversibility.  Following  their  models,  irreversibility  is  incorporated  in  the  present  work  by 
allowing  the  cohesive  surfaces  to  linearly  unload  to  and  reload  from  the  origin. 


3.2.2  Implementation  of  Cohesive  Zone  Models  in  VCFEM 

For  implementing  the  interfacial  constitutive  relations  in  VCFEM,  the  last  two  terms  in  equation 
4.21  are  replaced  with  traction-separation  relations  of  the  cohesive  zone  model.  For  Needleman’s 
models  this  takes  the  form 


-  /  r 


Un  +  AUn 


TTdun  ddo. 


-I  r 

Jut 


Uf  +  Aut 


T^dut  ddn 


=  /  [(^(Un -I- AUn,  Ut  -I- Alt()  -  C)(Un.  U()jd5n 

Jduum 


(3.16) 


where  the  potentials  (^(unjUt)  are  stated  in  equations  (3.9.3.11).  For  Geubelle’s  model  the  work  of 
separation  becomes 


J^dun  ddn 


11 


/  / 

(1-^)  / 
Jd 


ut  4-Atit 


dan 

—  ^(<)n  +  +  i^t  +  ddfi 


an^/an?  2  Smai 

c 


^  A'  +  (Sn  +  A5„)2  +  Jsi  +  {5t  +  ASt)^ 

Van^/ang  l  -  ^max 

-  i  ((<5„  +  A(5„)2  +  (5t  +  A<Jt)2)]  dan 


(3.17) 


where  is  a  step  function  defined  as 


'i  = 


0  for  S  ^  ^max 
1  for  6  >  ^max 


The  normal  and  tangential  separation  at  the  interface  may  be  expressed  in  terms  of  the  nodal 
displacements  using  the  interpolated  forms  of  equation  4.22  as: 

„„=<-<  =  {n'}’‘[L1{q'"  -  q‘} 
u,=  ur-uf  =  {t‘)’‘[L'|{q'"-q‘}  (3.18) 

Here  {n^^}  and  {t*^}  are  arrays  containing  direction  cosines  of  unit  outward  normal  and  unit  tangent 
vectors  to  the  interface.  The  resulting  weak  forms  of  the  element  kinematic  relations  and  traction 
reciprocity  conditions  are  nonlinear  due  to  the  nonlinear  traction-displacement  relation  in  the  cohe¬ 
sive  zone  models.  A  Newton-Raphson  method  iteration  method  is  consequently  employed  to  solve 
for  the  stresses  and  displacements.  In  the  event  of  unloading  in  an  increment  due  to  reduction  of 
interfacial  separation  (crack  closure),  an  unloading  algorithm  to  approach  the  origin,  is  activated. 


3.3  Numerical  Examples 

The  numerical  examples  are  divided  into  three  categories.  In  the  first  example,  VCFEM  results 
with  the  cohesive  zone  model  are  compared  with  results  in  Needleman  [22].  The  second  set  of 
examples  is  aimed  at  the  evaluation  of  cohesive  parameters  by  comparison  of  VCFEM  simulations 
with  experimental  results  and  further  validation.  In  the  final  set,  VCFEM  simulations  are  used 
to  make  predictions  of  damage  and  overall  mechanical  behavior  for  fiber  reinforced  composite 
microstructures  with  different  morphologies.  The  stress  functions  in  the  inclusion  phase  of  each 
Voronoi  Cell  element  is  generated  using  33  terms  (7-th  order  polynomial  stress  function,  i.e.  p+q  = 
2. .7)  for  the  polynomial  function  in  equation  (3.4).  The  matrix  stress  function  has  an  additional  36 
reciprocal  terms  due  to  the  reciprocal  terms  in  equation  3.5  (3  reciprocal  terms  for  each  polynomial 
exponent  from  2  to  4,  i.e.  i  =  p  +  q..p  +  q  +  2^p  +  q  €  [2.4]).  Displacement  fields  on  the 
•element  boundary  cind  on  the  matrix  and  inclusion  parts  of  the  interface  are  represented  using  linear 
interpolations  [L*]  and  [L*^]  in  equation  (4.22).  The  number  of  node-pairs  on  the  interface  is  set  to 
16.  This  is  consistent  with  the  requirements  of  stability  and  rank  sufficiency  of  the  assumed  stress 
hybrid  FEM,  as  discussed  in  [34,  21).  The  initial  locations  of  interfacial  nodes,  prior  to  debonding, 
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are  determined  by  an  adaptation  procedure  that  is  detailed  in  [21].  In  this  process,  nodes  are  added 
or  positioned  to  minimize  the  error  in  inter-element  and  interfacial  traction  continuity.  It  results 
in  optimization  of  the  virtual  work  with  respect  to  traction  discontinuity  across  boundaries.  The 
node-pairs  are  moved  in  the  direction  of  the  debond  once  debonding  sets  in,  to  accurately  capture 
the  crack  tip  stress  concentration. 


3.3.1  Comparison  with  results  of  Needleman  [22] 

In  [22],  the  cohesive  zone  model  of  equations  (3.11)  has  been  used  for  simulating  debonding  of 
an  elastic-plastic  matrix  from  rigid  spherical  reinforcements.  The  reinforcements  of  1.04%  volume 
fraction  are  arranged  in  a  cylindrical  periodic  array.  Consequently,  an  axisymmetric  analysis  of  a 
square  unit  cell  with  periodic  boundary  conditions  is  conducted.  In  Needleman’s  analysis  [22],  the 
matrix  is  modeled  as  rate  dependent  elastic-viscoplastic  material  with  isotropic  hardening,  with 
the  effective  plastic  strain  rate  e  being  expressed  in  terms  of  the  effective  stress  a  as 

^(e)  =  cTo{l  +  e/eo)'*  •  cq  =  —  (3.19) 

where  m  and  N  are  the  strain  rate  hardening  and  strain  hardening  exponents  and  ctq  is  a  reference 
strength.  The  matrix  material  properties  are  assumed  to  be  Em  =  SOOctq,  i^m  =  0.3  and  ctq  = 
400M Pa.  The  current  version  of  VCFEM  formulation  can  only  support  rate  independent  plastic¬ 
ity.  It  is  deemed  that  the  small  value  of  the  exponent  (m  =  0.01)  in  [22]  would  agree  (at  .least 
qualitatively)  with  the  VCFEM  results  using  rate  independent  plasticity  (m  =  0).  The  inclusion  is 
modeled  as  elastic  with  very  high  stiffness.  The  cohesive  zone  parameters  in  equations  (3.11)  are 
taken  as  tT^ai  =  3<to,  S/ro  =  0.01,  and  a  =  10.0.  The  stress  triaxiality  parameter  p,  is  set  to  0.5. 
The  VCFEM  model  consists  of  a  single  element  with  periodic  boundary  conditions  imposed  as 


on  X  =  0  Ux  =  Ux  =  €ccf>o> 

ily  —  Oy, 


on  X  =  Ro 

on  y  =  Rq  (3.20) 


Here,  lii  and  Uy,  and  Tj  and  ty  are  the  velocities  and  traction  rates,  in  the  radial  and  axial 
directions  respectively.  The  macroscopic  effective  stress  Ee  =  ISi  -  Syi  is  plotted  as  a  function 
of  the  axial  strain  Ca  =  /n(l  4-  Ux/bo)  in  figure  3.2.  The  onset  of  debonding  is  signaled  by  the 
sudden  drop  in  stress  level  and  its  arrest  is  indicated  by  the  resumption  of  monotonic  increase  of 
flow  stress.  The  results  of  VCFEM  simulation  agree  rather  well  with  those  in  [22].  The  slightly 
higher  stress  values  can  be  attributed  to  the  small  deformation  rate  independent  plasticity  with 
exponent  m  =  0  in  VCFEM  formulation. 


3.3.2  An  Experimental-Computational  Study 

Prior  to  predicting  the  interfacial  debonding  phenomena  for  multiple-fiber  reinforced  composites 
with  the  computational  model,  limited  comptirisons  with  experiments  are  conducted  with  model 
composite  systems  with  two  objectives.  The  first  is  to  ju-^certain  the  material  parameters  in  the 
cohesive  zone  model,  viz.  (Tmaxi  ttnd  oi  in  Needleman  s  models  [24]  and  (Xmaxi  ^maxi  ^’iid 
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u?  in  Geubelle’s  model  [15].  The  parameters  are  evaluated  by  solving  inverse  problems  in  which 
the  diflFerence  between  critical  experimental  observations  and  results  of  VCFEM  simulations  are 
minimized.  The  second  objective  is  to  validate  the  predictions  of  VCFEM  against  experimental 
observations  for  composites  with  different  microscopic  architectures. 

The  experimental  procedure 

The  debonding  experiments  are  conducted  with  model  single  and  multiple  fiber  specimens  in  the 
form  of  a  cruciform  as  shown  in  figure  3.3a;  The  cruciform  shape  has  been  developed  in  [12,  19,  28] 
to  avoid  the  development  of  stress  singularity  at  the  intersection  of  fiber-matrix  interface  and  free 
surface,  with  conventional  uniform  width  specimens.  Such  stress  singularities  promote  interfacial 
separation  near  the  free  surface  and  results  in  invalid  interface  strength  data.  The  most  significant 
advantage  of  the  cruciform  geometry  is  that  it  forces  debond  failure  to  initiate  in  the  central  region 
of  the  specimen. 

The  model  composite  specimens  are  fabricated  by  casting  in  a  cruciform-shaped  silicone  rubber 
mold  to  specimen  dimensions  shown  in  figure  3.3  and  table  1.  The  circular  fillets  at  the  cross 
junctions  reduce  stress  concentration  in  the  matrix  at  these  locations.  The  reinforcing  fibers  in 
the  composite  system  are  stainless  steel  filaments.  The  fibers  of  predetermined  length  are  posi¬ 
tioned  in  the  mold  and  the  matrix  material  is  cast  around  them.  Larger  diameters  are  helpful 
in  detecting  the  debond  initiation  process  and  are  chosen  for  the  single  .fiber  cruciforms.  The 
matrix  is  an  epoxy  resin  (Epon  828  manufactured  by  Shell  Chemical  Co.),  that  is  cured  with  a 
polyetheramine  (Jeffamine  D-230  manufactured  by  Texaco  Inc.)  for  3  days  at  ambient  temperature. 

Room-temperature  curing  reduces  residual  thermal  stresses  in  the  matrix,  due  to  mismatch  in 
the  fiber  and  matrix  thermal  properties,  to  a  minimum.  The  epoxy  matrix  is  transparent  and  al¬ 
lows  visualization  of  the  debonding  process  at  the  fiber-matrix  interface.  Preliminary  experiments 
with  no  interfaeial  coating,  indicate  that  the  bond  strength  is  high  and  interface  failure  progresses 
rapidly  to  cause  immediate  catastrophic  rupture.  Consequently,  the  steel  filaments  are  cleaned  and 
polished  with  acetone  and  coated  with  a  very  thin  film  of  freekote  (<  0.1  y.m)  prior  to  casting 
the  specimens.  The  freekote  imparts  a  weak  strength  to  the  steel-epoxy  interface.  This  allows  a 
somewhat  stable  growth  of  the  debond  crack,  thereby  permitting  determination  of  the  parameters 
in  interfacial  cohesive  zone  models. 

In  the  experimental  set  up,  three  strain  gages  (A,  B  and  C)  are  affixed  on  faces  of  the  specimen 
as  shown  in  figure  3.3a.  Two  gages  (A  and  B)  are  locattnl  in  the  central  portion  of  the  cruciform  in 
vicinity  of  potentitil  debond  sites.  Thus  their  readings  are  .ussnined  to  correspond  to  the  debonding 
strains  in  the  specimen.  The  third  gage  C  is  mounted  on  the  limb,  away  from  the  fibers  and 
represents  the  fax-field  strain.  Furthermore,  to  prevent  spe(  imeti  failure  in  the  gr^p  region,  fiber¬ 
glass/epoxy  end  tabs  are  adhesive  bonded  on  the  upright  portion  of  the  specimen.  Four  different 
microstructural  architectures  are  considered  for  the  experiments.  They  are: 

•  Specimen  #  1  containing  a  single  circular  fiber. 
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•  Specimen  #  2  containing  a  single  elliptical  Eber  with  major  axis  along  the  loading  direction 

(Omaj  —  0°), 

•  Specimen  #  3  containing  a  single  elliptical  fiber  with  major  axis  perpendicular  to  the  loading 
direction  {amaj  =  90°), 

•  Specimen  #  4  containing  5  identical  circular  fibers. 

These  are  schematically  illustrated  in  figure  3.3  and  the  geometric  details  on  location,  size  etc. 
are  given  in  table  1.  The  area  fraction  in  the  table  are  for  the  cross-section  at  the  outer  edge. 
The'elastic  material  properties  for  the  reinforcing  steel  fibers  and  epoxy  matrix  are  experimentally 
determined  as: 

Young^s  modulus  Esteel  =  210  GPa,  Poisson's  ratio  Usteei  =  0.3;  Young’s  modulus  Eepoxy  =  4.6 
GPa,  Poisson’s  ratio  ^epoxy  — 

The  cohesive  properties  of  the  debonding  material  are  evaluated  by  solving  the  inverse  problem 
mentioned  in  section  3.3.2. 

The  model  specimens  are  loaded  in  tension  on  a  servo- hydraulic  testing  machine.  The  readings 
from  the  strain  gages  are  recorded  continuously  and  the  onset  of  fiber-matrix  debonding  is  inter¬ 
preted  from  a  sharp  change  in  the  slope  of  the  stress-strain  curve  based  on  gage  readings  at  the 
cruciform  center.  Acoustic  emission  sensors  are  also  employed  to  confirm  the  onset  of  debonding. 
In  figure  3.5,  the  abrupt  change  in  slope  at  B  corresponds  to  debond  initiation,  and  the  relatively 
flat  portion  BC  corresponds  to  the  strain  jump.  Subsequent  loading  (CD)  proceed  with  a  lower 
stress-strain  slope,  due  to  a  reduced  load  caring  capability  of  the  partially  debonded  fiber.  Un¬ 
loading  along  DA  and  subsequent  reloading  along  ACD  indicate  no  further  change  in  slope,  i.e. 
no  additional  debonding.  Matrix  cracks  always  initiate  at  the  fiber-matrix  debond  site  and  grow 
rapidly  to  cause  specimen  failure.  A  few  specimens  are  not  loaded  to  fracture,  to  allow  observation 
of  the  partially  debonded  interface. 

Following  interface  failure,  some  of  the  cruciform  specimens  are  sectioned  along  the  center, 
parallel  to  the  loading  direction.  A  drop  of  a  fluorescent  dye  penetrant  is  positioned  above  the 
sectioned  fiber  and  vacuum  infiltration  is  used  to  force  the  penetrant  into  the  debonded  interface. 
The  sectioned  face  is  then  polished  to  remove  traces  of  the  dye  from  the  original  drop.  Thereby,  only 
the  dye  that  has  penetrated  along  the  interface  crack  will  remain.  The  fiber  is  then  viewed  under 
an  ultraviolet  light,  which  cause  the  regions  containing  the  dye  to  remain  bright  in  an  otherwise 
dark  background.  Figure  3.4a  shows  the  debond  under  oblique  incidence,  where  the  fiber  surface 
is  visible  through  the  transparent  matrix.  The  bright  regions  correspond  to  locations  of  debond 
and  they  are  concentrated  on  the  fiber  surface  in  the  direction  of  the  loading  axis.  Also,  note  that 
they  are  generally  absent  in  the  wing  region,  which  experience  significantly  lower  loads  than  the 
cruciform  center.  Figure  3.4b  shows  the  cross  section  of  the  fiber  at  a  higher  magnification,  and 
the  debonded  region  is  the  thin  bright  strip  along  the  fiber  periphery.  The  ends  of  the  debond  are 
highlighted  by  the  arrows,  and  the  loading  direction  is  horizontal  in  the  figure.  Figure  3.4b  also 
shows  that  the  total  angle  of  debonding  was  approximately  85°. 
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Evaluation  of  a  Conversion  Factor  for  Stress  3-D  to  2D  Stresses 

The  experimental  results  are  compared  with  results  of  simulations  by  a  2-D  (plane  strain)  Voronoi 
cell  finite  element  model  for  validation  and  damage  prediction.  While  the  fiber- reinforced  cruciform 
specimen  represents  a  3-D  problem  in  actuality,  the  stresses  and  displacements  at  the  center  of  the 
cruciform  can  be  well  approximated  by  two  2-D  solutions  for  a  reasonably  slender  geometry.  Prior 
to  comparison  of  the  results,  a  conversion  factor  is  established  for  predicting  the  3-D  stress  state 
from  2-D  simulations.  A  full  3-D  analysis  of  the  cruciform  specimen  and  a  plane  strain  analysis 
of  the  section  containing  a  single  circular  fiber  are  conducted  using  the  commercial  finite  element 
package  ANSYS.  Perfect  interface  with  no  debonding  is  assumed  in  these  analyses.  Figu-e  3.6 
shows  comparison  of  the  radial  stress  Orr  and  tangential  stress  are  along  the  interface  for  the  two 
analyses.  A  stress  concentration  factor  is  defined  as  the  ratio  of  the  maximum  radial  stress  along 
the  tensile  axis  at  the  pole  of  the  fiber-matrix  interface  (angle=0  in  the  figure)  to  the  far-field 
stress.  This  factor  is  computed  to  be  k2D  =  1-349  for  2D  analysis  and  k^o  -  1-139  for  3-D.  By 
multiplying  the  plane  strain  results  with  a  factor  k  =  ^  =  0.844  enables  a  very  good  match  with 
the  3-D  solution,  as  shown  in  the  stress  plots  of  figure  3.6.  Consequently,  the  2-D  stresses  from 
VCFEM  analysis  are  multiplied  by  the  factor  k  =  0.844  for  comparison  with  3-D  experimental 
results.  The  2-D  sections  analyzed,  corresponds  to  rectangular  vertical  sections  through  the  center 
of  the  cruciform  specimen,  shown  in  figure  3.3a.  The  sections  have  dimensions  of  33  mm  x  ye  with 
the  fiber  cross-sections  in  the  central  region,  and  where  Ve  is  given  in  table  1.  A  total  of  (n  -f-  2) 
Voronoi  cell  elements,  where  n  correspond  to  the  fibers  and  two  additional  elements  for  the  limb 
portions,  are  used  to  perform  this  2-D  analysis. 

Evaluation  of  cohesive  zone  parameters 

Cohesive  zone  parameters  for  the  models  proposed  by  Needleman  [22,  23]  and  Geubelle  et.  al. 
[9.  15]  are  calibrated  by  solving  an  inverse  problem  using  the  VCFE  model  and  comparing  with 
experimental  results.  The  computational  model  analyzes  the  2-D  cross-sections  shown  in  figure 
3.3b  under  plane  strain  conditions  and  loaded  axially.  In  Geubelle’s  model,  the  parameters  to  be 
determined  are  6max,  cTmax^  u?  and  u^,  while  in  Needleman's  model  the  parameters  are  CTmai.  and 
Q.  Their  evaluation  for  a  given  bonding  material  (freekote  in  this  problem)  follows  a  multi  variable 
optimization  process,  in  which,  four  key  variables  are  chosen  for  minimizing  the  difference  between 
experimental  and  computational  results.  These  effective  values  are  marked  in  figure  3.5  and  are; 
(i)  the  strain  at  which  debonding  initiates  (5/),  (ii)  the  strain  (5//)  at  which  debonding  arrests, 
(iii)  the  stress  (5/;/)  at  which  debonding  arrests,  and  (iv)  post  debonding  slope  of  the  stress-strain 
plot  (Siv).  The  major  steps  in  this  evaluation  process  are; 

1.  For  each  of  the  models  in  equations  3.11,  3.13  and  3.13.  the  overall  stress-strain  response  with 
debonding  behavior  is  first  determined  (similar  to  figure  3.5)  for  a  range  of  candidate  values  of 
cohesive  parameters.  These  parameters  cire  chosen  from  an  estimated  range  c£possible  values. 
For  each  set  of  parameters,  e.g.  {6maxy  (^max^  ^n)  o'"  {<^rnax^  <5’  and  J*),  an  independent 
VCFEM  simulation  of  the  microstructure  with  debonding  in  simple  tension,  is  performed. 

2.  The  key  variables  Si{i  —  I...IV)  are  represented  as  polynomial  functions  of  the  cohesive  zone 


16 


parameters.  This  produces  continuous  functions  for  gradient  based  optimization  methods, 
necessary  in  parameter  evaluation.  For  Geubelle’s  model  [15,  9],  these  are  written  as 

Si  =  (Ojo  +  +  at2<^mai +  •••)  + (^lO  +  +  •••) 

+  {cio  +  CiiUi  +  Ci2Uf^ +  ...)  + [dio  +  diiu^  +  di2U^^  +  ...)  (3.21) 

while  for  Needleman’s  model  [22,  23],  these  are 

Si  =  (oio  +  anO’mox  +  fli2<7max  +  ••)  +  (^*0  +  +  i)t2<^*^  +  ••)  +  (ctO  +  cua  +  caoP'  +  .CB-22) 


3.  To  evaluate  the  coefficients  of  the  polynomial  expansions  i.e.  Ujo,  an , ...,  bio,  bn , ...,  Cjo,  Cji , ..., 
dio.dti,...,  approximately  150  VCFEM  simulations  are  conducted  with  different  parameter 
sets.  For  each  simulation  with  a  set  of  values  e.g.  (Smax,  amaxi  u^)  or  (amoxi  <5*  S'lid 
(5*),  the  computed  values  of  Si,i  =  I. .IV  are  recorded  from  the  stress-strain  plots.  The 
minimum  number  of  simulations  correspond  to  the  total  number  of  unknown  coefficients. 
However,  a  higher  number  is  performed  in  this  study  to  accommodate  a  larger  range  of 
cohesive  parameters. 

4.  The  coefficients  aio,Oti,.-M  bio,bn,...,  Cio,cn,...,  diQ,dn,...  are  evaluated  by  a  least  square 
based  minimization  process  using  MATLAB.  This  enables  explicit  construction  of  the  func¬ 
tions  S,,i  =  I. .IV. 

5.  The  parameters  of  cohesive  zone  models  in  Geubelle  and  Needleman  are  finally  obtained  by 
minimizing  the  difference  between  experimental  and  simulated  values  of  the  key  vairiables  Si 
in  the  stress-strain  plots  of  figure  3.5  or  3.10.  A  multi-parameter  optimization  problem  is 
solved  as  : 


Minimize 


±wi 

t=l 


i(s? 


simulation 


W.T.t.  (Tmaxi  ^max’t  ^71 

OR 


^tnax  i  ^  OLTtdt  Ol 


gexperiment^ 


2 


(3.23) 


The  weights  w,  may  be  suitably  chosen  to  impart  preferential  importance  to  the  key  variables. 

Equal  weights  au'e  chosen  for  this  analysis. 

Two  sets  of  resulting  values  of  the  cohesive  zone  parameters  are  presented  in  table  2.  The  first 
set  correspond  to  values  that  are  evaluated  by  comparison  with  the  results  of  experiments  with 
specimen  1.  When  these  parameters  are  used  to  simulate  the  other  specimens  2,  3  and  4  ,  it  is  seen 
that  the  maximum  difference  between  experiments  and  simulations  occur  for  specknen  3.  A  second 
set  of  parameters  are  consequently  derived  using  the  key  variables  Si  for  both  specimens  1  and  3, 
shown  in  figures  3.5  and  3.10.  A  small  difference  is  noted  in  table  2,  between  the  values  of  the  two 
sets.  The  results  shown  in  the  following  sections  are  with  the  second  set  of  cohesive  parameters. 
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Comparison  of  VCFEM  and  experimental  results 

Results  of  VCFEM  simulations  for  the  four  different  microstructures  of  specimens  1,  2,  3  and  4  are 
compared  with  experiments,  macroscopically  in  figures  3.5,  3.10  and  3.15  and  microscopically  in 
figure  3.4.  The  abscissa  in  figures  3.5,  3.10  and  3.15,  is  the  gage  strain  recorded  by  0.8  mm  strain 
gages  mounted  on  the  specimen  surface  (A  in  figure  3.3a).  The  corresponding  strains  in  VCFEM 
analyses  are  calculated  from  the  change  in  length  of  a  0.8  mm  segment,  located  on  the  2D  specimen 
boundary,  closest  to  the  fiber.  In  table  3,  the  simulated  gage  strains  and  the  area-averaged  strains 
are  tabulated.  There  is  a  small  difference  between  the  two  strains.  This  may  be  attributed  to  the 
fact  that  the  edge  along  which  the  gage  strain  is  measured  is  sufficiently  close  to  the  fibers  for 
local  effects  of  the  interface  to  prevail.  The  ordinate  in  the  plots  of  figures  3.5,  3.10  and  3.15  are 
the  area-averaged  macroscopic  stresses  in  the  direction  of  loading.  The  stress-strain  plots  for  all 
three  decohesion  models  (Geubelle,  Needleman-polynomial  and  Needleman-exponential)  in  figure 
3.5  show  good  agreement  with  the  experimental  stress-strain  behavior.  The  onset  of  debond,  sig¬ 
naled  by  a  sharp  reduction  of  slope,  is  also  predicted  quite  well.  However,  the  arrest  of  debond,  as 
indicated  by  an  increase  in  the  post-decohesion  stiffness,  is  more  gradual  for  the  simulations  than 
for  experiments. 

The  experimentally  determined  debonding  angle  by  dye  injection  is  evaluated  in  figure  3.4  to 
be  approximately  85°.  The  corresponding  simulated  value  is  found  to  be  91°  (  see  table  3).  This 
corresponds  to  a  6.6%  deviation  from  experiments.  Table  3  also  presents  the  overall  stiffness  values 
before  and  after  debonding  and  the  total  debonding  angle  along  the  interface,  using  Geubelle’s 
model  [15].  The  undamaged  overall  stiffness  K^d  is  the  slope  of  the  area  averaged  tensile  stress- 
strain  plots  in  the  loading  direction,  for  the  rectangular  section  through  the  cruciform  center.  The 
highest  K^d  is  observed  for  specimen  2  with  a^aj  =  0°,  while  the  lowest  is  that  for  specimen  4 
with  a  lower  volume  fraction.  The  maximum  drop  in  stiffness  occurs  with  specimen  3  {ocmaj  ~  90°) 
since  this  interface  undergoes  the  largest  debonding  angle  (0^=134°).  The  minimum  drop  is  for  the 
multiple  fiber  specimen  even  with  considerable  interfacial  debonding  (^<i=87°),  mainly  due  to  the 
small  area  fraction.  The  increase  in  debonding  angle,  as  functions  of  the  averaged  strain,  are  plot¬ 
ted  in  figures  3.7  and  3.16.  For  specimen  1,  the  response  of  the  three  interfacial  decohesion  models 
are  extremely  close.  The  debonding  angle  increases  rapidly  after  initiation  with  little  additional 
strain,  and  then  tapers  off  asymptotically  to  a  stable  value.  The  initiation  strains  for  specimens  1 
and  3  are  quite  close,  while  that  for  the  specimen  2  is  significantly  lower.  Also  the  smallest  stable 
debond  angle  occurs  with  the  specimen  2  and  the  largest  for  specimen  3  with  specimen  1  in  the 
middle.  The  stable  debonding  angle  for  both  specimens  I  and  4  with  circular  fibers  are  quite  close. 

To  understand  the  effects  of  stress  distribution  on  tlu'  ilebonding  process,  the  normal  and  tan¬ 
gential  interfacial  stresses  CTnn  and  Cnt  are  plotted  from  o  =  0"  (the  loading  direction)  to  4>  =  180°, 
in  figures  3.8,  3.9,3.11,3.12,  3.13  and  3.14.  The  peak  stress  at  the  tip  of  the  crack  or  debond 
is  obtained  from  the  spring  tractions  at  the  node-pair,  just  ahead  of  the  debonefed  nodes  where 
the  condition  5  =  1  in  equation  (3.14)  is  satisfied.  The  normal  stress  plots  are  symmetric  about 
the  angular  position  <p  —  90°,  while  the  tangential  stresses  are  antisymmetric  about  this  angle. 
Prior  to  debonding,  the  normal  stress  is  mriximum  at  o  =  0°.  while  the  shear  stress  peaks  out  at 
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4)  =  45“.  Once  debonding  initiates,  the  peaks  occurs  at  the  crack  tip  and  subsides  monotonically  till 
it  reaches  the  pole  at  0  =  90",  where  is  mildly  compressive.  With  progressive  debonding,  the 
peak  tensile  stress  (7„„  initially  increases  in  magnitude  but  subsequently  decreases.  This  behavior 
may  be  explained  as  a  consequence  of  two  competing  phenomena,  viz.  an  increase  in  average  stress 
due  to  increasing  debond  length  and  a  decrease  in  the  normal  component  of  stress  with  increasing 
angular  orientation,  especially  as  90°.  The  compressive  region  also  increases  with  increasing 
decohesion.  Similar  observations  have  also  made  in  [1].  For  the  tangential  stress  ant,  the  maximum 
value  at  the  crack  tip  is  also  found  to  fiist  increase  slightly  and  then  decrease  with  progressive 
debonding.  From  these  plots,  it  is  easy  to  explain  the  reasons  for  the  nucleation  and  progress  in 
decohesion,  as  depicted  in  figures  3.7  and  3.16. 

Through  a  combination  of  experimental  and  computational  simulations,  these  examples  estab¬ 
lish  considerable  confidence  in  the  model  for  understanding  the  behavior  of  the  interfacial  debonding 
phenomena. 

3.3.3  Multiple  Fibers  in  the  Microstructure 

In  this  final  set  of  examples,  computational  analyses  are  conducted  to  understand  the  effect  of 
morphology  and  applied  boundary  conditions  on  the  initiation  and  progress  of  debonding  induced 
damage  in  microstructures  containing  multiple  fibers.  The  microstructure  in  these  examples  contain 
100  circular  fibers  in  a  square  region.  To  investigate  the  role  of  fiber  interaction  due  to  different 
spatial  dispersions,  two  distributions  viz.  a  uniform  distribution  and  a  random  distribution  are 
considered.  For  comprehending  the  role  of  boundary  conditions  on  the  damage  evolution  process, 
two  in-plane  conditions  viz.  (a)  periodic  conditions  and  (b)  displacement  (non-periodic)  boundary 
conditions,  are  specified.  The  periodicity  boundary  conditions  are  imposed  by  requiring  edges  to 
remain  straight  and  parallel  to  the  original  direction  throughout  deformation.  This  may  be  achieved 
through  the  following  conditions. 

Uj  =  0  {on  X  =  0)  ,  Uy  =  0  {on  y  =  0)  ,  Uj  =  Uap  {on  x  =  Lx)  ,  Uy  =  Dy  {on  y  =  Ly) 

Ty  =  0  (oni  =  0  and  X  =  Lx)  ,  Ti  =  0  (on  y  =  0  and  y  =  Ly)  (3.24) 

where  Uap  is  a  monotonically  increasing  applied  displacement  and  D*  is  determined  from  the  av¬ 
erage  force  condition  J^Txdx  =  0  on  y  =  Ly.  For  the  displacement  boundary  conditions,  y  =  Ly 
is  a  traction  free-edge.  The  elastic  material  properties  for  the  reinforcing  steel  fibers  and  epoxy 
matrix  are  the  same  as  in  the  previous  section,  i.e.  'loung  s  modulus  E steel  =  210  GPa,  Poisson  s 
ratio  Ojtee;  =  0.3;  Young’s  modulus  Eepoxy  =  4.6  GPa.  Poisson  s  ratio  t^epoiy  =  0-4.  The  interface 
is  represented  by  the  cohesive  zone  model  of  Geubelle  ct.  al.  45]  with  the  same  parameters,  i.e. 
<5mai  =  0.9347,  Omcx  =  0.012570  GPa,  ul  =  0.00007217  mm  and  u°  =  0.00006395  mm.  The 
plane  strain  Voronoi  cell  finite  element  model  is  used  in  this  analysis,  with  a  tessellated  mesh  of 
100  Voronoi  elements  each  containing  a  fiber.  For  the  uniform  distribution,  the  tessellation  process 
yields  square  elements,  while  for  the  random  distribution  the  cells  have  variable  number  of  edges. 
The  two  meshes  with  boundary  conditions  are  shown  in  figures  3.17  a  and  b. 
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The  macroscopic  or  area-averaged  stress-strain  response  for  the  two  architectures  are  plotted 
for  the  two  boundary  conditions  in  figure  3.18a.  Prior  to  debonding,  the  stiffness  of  the  random 
microstructure  is  marginally  higher  than  that  for  the  uniform  one.  For  each  microstructure,  the 
behavior  with  periodicity  is  stiffer  that  that  without.  This  is  due  to  the  additional  displacement 
constraint  imposed  by  the  periodicity  condition,  which  raises  stresses  in  the  microstructure.  Appre¬ 
ciable  debonding  in  fibers  is  signaled  by  an  abrupt  change  in  slope  of  the  stress-strain  curve.  This  is 
noticed  at  a  significantly  smaller  strain  for  the  random  microstructure.  However,  the  ultimate  loss 
of  stiffness  is  higher  for  the  uniform  microstructure.  The  number  cf  debonded  fibers  as  a  function 
of  increasing  strain  is  plotted  in  figure  3.18b.  The  rate  of  damage  growth  is  more  gradual  for  the 
random  microstructure.  This  rate  is  initially  low,  but  increases  considerably  with  straining  in  the 
intermediate  region  and  finally  tapers  off  at  higher  values  of  strain.  This  is  because,  initial  stages  of 
debonding  are  influenced  by  the  inter-fiber  spacing.  Interfaces  of  fibers  that  are  in  close  proximity 
are  damaged  early.  The  subsequent  stages  of  damage  also  depend  on  the  overall  strain  state.  In  the 
uniform  model,  the  initiation  takes  place  later  but  the  growth  of  damage  by  interfacial  decohesion 
occurs  rapidly.  This  is  particularly  true  for  the  periodic  boundary  condition.  With  the  displacement 
boundary  conditions,  the  growth  is  slow  in  the  beginning  but  increases  substantially  with  straining. 

A  clear  delineation  of  the  debonding  process  is  obtained  from  the  contour  plots  of  the  tensile 
stress  axi  in  figures  3.19  and  3.20.  The  figures  show  the  debonds  on  the  interfaces,  as  well  as 
locations  of  stress  concentrations.  For  the  uniform  microstructure  with  displacement  boundary 
conditions  (non-periodic),  debonding  starts  at  the  fiber  in  the  upper  right  corner  and  propagates 
down  along  the  diagonal  to  create  a  dominant  damage  path.  The  reason  for  the  damage  propa¬ 
gation  along  a  distinct  path  is  explained  from  the  radial  stress  Orr  plots  along  the  interfaces,  in 
figures  3.21.  The  fiber  numbers  1,  2,  3,  and  4  are  delineated  in  figure  3.17a.  At  a  macroscopic 
strain  value  of  e  =  1.31  x  10"'*,  just  before  initiation,  the  radial  stress  at  the  interface  of  the  corner 
fiber  #2  is  higher  at  the  poles  along  to  the  loading  axis  (i.e.  d  =  0°,  180°).  This  causes  the  fiber 
#  2  to  debond  first.  Following  this,  the  radial  stress  in  fiber  #  3  exceeds  that  of  others  at  the 
poles  and  hence  debonds.  This  pattern  is  observed  for  the  entire  load  cycle.  Finally,  all  fibers 
debond  at  higher  values  of  the  applied  strain.  In  contrast,  all  fibers  debond  simultaneously  for 
the  periodic  boundary  conditions,  with  no  apparent  dominant  path.  A  single  unit  cell  with  these 
periodic  boundary  conditions  also  produce  identical  results.  It  is  clear  from  this  example,  that 
even  though  the  microstructure  may  have  geometric  periodicity,  the  evolution  of  damage  to  cause 
failure  is  likely  to  follow  a  non-periodic  and  dominant  local  path.  This  makes  the  application  of 
the  conventional  unit  cell  analysis  to  problems  with  evolving  damage,  inappropriate. 

For  the  random  microstructure,  the  difference  in  debond  initiation  and  growth  behavior  with 
different  boundary  conditions  is  not  as  prominent.  In  contrast,  the  effect  of  local  morphology  is 
much  more  dominant.  This  is  observed  in  figures  3.18b  and  3.20.  Debonding  initiates  in  fibers  that 
are  in  close  proximity  with  other  fibers  and  evolves  randomly  with  increased  strainmg.  The  relation 
between  damage  evolution  and  the  local  dispersion  for  the  random  microstructure,  is  better  under¬ 
stood  from  the  histograms  of  figures  3.22  and  3.23.  In  these  histograms,  the  relative  spacing  and 
proximity  among  fibers  is  characterized  by  two  spatial  functions  viz.  (a)  local  area  fraction  (LAF) 
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and  (b)  the  nearest  neighbor  distance  (NND).  The  local  area  fraction  (LAP)  is  measured  as  the 
ratio  of  the  fiber  area  to  the  area  of  the  Voronoi  cell  containing  it.  The  nearest  neighbor  distance 
(NND)  is  measured  as  the  smallest  surface  to  surface  distance  between  neighbors  that  share  com¬ 
mon  Voronoi  cell  edges.  The  ordinate  of  the  histograms  represents  the  total  number,  as  well  as  the 
number  of  debonded  fibers  as  functions  of  the  characterization  functions  LAP  and  NND.  Different 
shades  are  used  to  delineate  additionally  debonded  fibers  with  increasing  strain.  The  LAP  function 
has  a  bell  shaped  distribution  whereas  the  NND  function  is  clearly  monotonically  decreasing  for 
the  total  number  (debonded+intact)  fibers.  Debonding  initiates  at  the  higher  values  of  LAP  and 
lower  values  of  NND.  There  is  a  monotonic  decrease  of  debonded  fibers  with  decreasing  LAP  and 
increasing  NND.  The  monotonicity  is  particularly  obvious  for  the  latter  function.  The  low  NND 
represents  closely  packed  fibers  which  are  more  likely  to  debond  early,  in  the  loading  cycle.  Similar 
responses  are  obtained  for  both  periodic  and  displacement  boundary  conditions.  The  statistics  of 
debonding  are  also  tabulated  in  table  4. 

The  effect  of  volume  or  area  fractions  are  also  investigated  through  a  second  area  fraction  of 
Vf  =  10%.  Both  constant  and  randomly  varying  sizes  of  randomly  dispersed  fibers  are  considered 
for  the  10%  microstructure.  The  macroscopic  responses  for  these  models  are  plotted  in  the  figures 
3.18.  As  long  as  the  volume  fraction  is  the  same,  the  overall  stiffness  is  not  very  sensitive  to  the 
size  distribution,  i.e.  it  is  quite  close  for  both  the  constant  and  the  randomly  varying  diameters. 
However,  the  dependence  of  damage  on  the  size  distribution  is  considerable  and  occurs  earlier  for 
the  variable  size.  The  onset  of  debonding  is  considerably  delayed  for  the  lower  volume  fraction, 
due  to  increased  average  spacing  between  fibers  which  result  in  lower  stress  concentration.  In 
summary,  the  process  of  initiation  and  growth  of  damage  by  interfacial  debonding  depends  on 
several  competing  factors  such  as  morphology,  fiber  interaction,  load  and  boundary  conditions.  It 
is  important  to  include  large  portions  of  the  microstructure  to  capture  the  different  effects,  which 
is  a  major  shortcoming  of  unit  cell  analyses. 


3.4  Concluding  Remarks 

The  initiation  and  growth  of  damage  by  interfacial  decohesion  in  multiple-fiber  polymer  matrbc 
composites  is  analyzed  with  in-plane  loading  in  this  work.  The  fiber  and  matrix  phases  are  mod¬ 
eled  with  elastic  properties.  The  progress  of  interfacial  debonding  with  quasi-static  loading  is 
modeled  by  cohesive  zone  constitutive  relations  in  terms  of  normal  and  tangential  tractions  and 
interfacial  separation.  In  these  relations,  the  traction  increases  with  sepairation,  reaches  a  maximum 
and  subsequently  subsides  to  zero  traction,  signaling  debonding.  The  stress  and  damage  analyses  is 
conducted  with  the  Voronoi  cell  finite  element  model,  that  has  been  established  as  an  effective  tool 
for  modeling  of  non-uniform  microstructures.  VCFEM  has  been  developed  (see  [34,  21])  to  yield 
high  computational  efificiency  with  good  accuracy  for  modeling  large  heterogeneousjnicrostructures. 

The  VCFE  model  for  interfacial  debonding  has  been  compared  with  a  variety  of  established 
results  in  literature,  the  details  of  which  are  reported  in  the  M.S.  thesis  of  Ling  [16].  Only  a  single 
satisfactory  comparison  of  macroscopic  results  with  those  in  [22]  are  presented  in  this  work.  A  com- 
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bined  experimental-computational  study  is  conducted  for  specimens  of  polymer  matrix  composites 
in  this  work  An  inverse  methodology  is  developed  for  calibrating  the  cohesive  model  parameters 
from  simple  experimental  results.  It  is  based  on  minimizing  the  difference  in  the  macroscopic  stress- 
strain  plots  with  damage,  obtained  by  VCFE  simulations  and  experiments.  For  the  experiments, 
single  and  multiple  fiber  cruciforms  specimens  of  circular  and  elliptical  shapes,  are  fabricated  to 
avoid  stress  concentrations  at  the  edges.  Results  are  considered  from  more  than  one  experimental 
specimens,  to  account  for  the  variability  between  experiments  in  the  calibration  process.  VCFE 
debonding  analyses  are  subsequently  conducted  to  simulate  experiments  with  other  specimens.  The 
microscopic  debonding  angle  in  experiments  is  estimated  by  a  dye  penetration  technique.  Good 
agreement  is  obtained  between  experiments  and  the  simulations,  both  with  respect  to  macroscopic 
(averaged  stress-strain  behavior)  and  microscopic  (debond  angle)  observations. 

In  numerical  simulations  of  multiple  fibers  microstructures,  different  morphologies  and  bound¬ 
ary  conditions  are  analyzed  to  understand  their  influence  on  the  decohesion  process.  Different 
architectures  include  uniform  and  random  dispersions,  different  volume  fractions  and  different  size 
distributions.  The  boundary  conditions  include  imposed  periodicity  and  prescribed  displacement 
conditions.  For  the  uniform  microstructure,  the  path  of  growing  damage  is  found  to  be  very  sen¬ 
sitive  to  the  boundary  conditions.  Even  with  periodic  geometric  features,  a  distinct  non-periodic 
and  dominant  damage  path  evolves  with  increasing  strain,  for  non-periodic  boundary  conditions. 
This  effect  is  significant  at  higher  volume  fractions.  Periodic  damage  is  only  observed  with  periodic 
conditions.  An  important  conclusion  that  can  be  derived  from  this  example  is  that  unit  cells  with 
periodic  boundary  conditions  are  inadequate  for  properly  predicting  the  expected  local  nature  of 
damage  growth.  For  random  microstructures,  the  debond  induced  damage  is  found  to  be  more 
sensitive  to  the  local  morphology  e.g.  inter-fiber  spacing  than  to  the  boundary  conditions.  The 
damage  growth  process  is  observed  to  take  place  gradually  over  a  longer  range  of  increasing  strain 
for  the  random  microstructures. 

This  study  reveals  the  significance  of  analyzing  large  regions  of  the  microstructure  and  proves 
the  effectiveness  of  the  VCFEM  analysis  for  the  same.  The  Voronoi  cells  also  play  an  important  role 
in  developing  geometric  descriptors  for  quantitative  characterization  (e.g.  LAF,  NND)  since  they 
represent  regions  of  immediate  influence  for  each  fiber  and  also  delineate  neighbors.  It  provides 
the  essential  link  between  the  microstructural  features  and  response,  that  is  important  in  damage 
analysis. 
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- 

4 

5 

Circular 

4.31,2.23 

7.65.5.69 

0.355.0.355 

1.14 

- 

4 

5 

Circular 

4.22,3.47 

7.65,5.69 

0.355,0.355 

1.14 

_ 

4 

5 

Circular 

2.94.3.37 

7.65.5.69 

0.355,0.355 

1.14  ■ 

- 

4 

5 

Circular 

3.68,2.78 

7.65,5.69 

0.355,0.355 

1.14 

- 

Table  3.1:  Statistics  of  geometrical  parameters  in  the  experimental  specimens:  (a)  specimen  num¬ 
ber.  (b)  number  of  fibers,  (c)  fiber  cross-section,  (d)  fiber  centroidal  coordinates  (xc,t/c),  (e)  size 
of  the  cross-section  shown  in  figure  3.1b  (ze,ye),  (f)  fiber  major  and  minor  axes  (a,  6),  (g)  cross- 
sectional  area  fraction  (A/)  and  (h)  major  axis  angle  with  the  loading  direction  (omaj)- 


Figure  3.1:  (a)  A  mesh  of  Voronoi  cell  elements  generated  by  tessellation  of  the  heterogeneous 
microstructural  domain  (b)  A  typical  Voronoi  cell  element  with  delineation  of  the  phases  and 
interface  region. 


26 


Model 

Specimen  # 

^max 

^max  1 

Geubelle 

1  only 

0.9412 

0.012846  GPa 

0.00007392  mm 

0.00006407  mm 

Geubelle 

1  k  3 

0.9347 

0.012570  GPa 

0.00007217  mm 

0.00006395  mm 

Model 

Specimen  # 

^max 

(5* 

a 

;  N'eedleman(poly) 

1  only 

0.014613  GPa 

0.00004837  mm 

•  Needleman(poly) 

1  k  3 

0.014428  GPa 

0.00004889  mm 

IIIIBQIQIIIIIIII 

i  Xeedleman(expo) 

1  only 

0.014980  GPa 

0.00005019  mm 

I  Needleman(expo) 

1  k  3 

0.014831  GPa 

0.00005007  mm 

0.6887 

Table  3.2:  Cohesive  zone  parameters  for  the  three  models  evaluated  from  the  specimen  1  only  and 
specimens  1  and  3  combined. 


Specimen  # 

/Cud(MPa) 

A:d(MPa) 

AK% 

Averaged  ef„% 

Gage 

^3 

1 

4974 

4265 

14.25 

0.140 

0.141 

91 

2 

5249 

4621 

11.96 

0.066 . 

0.063 

46 

3 

5163 

3581 

30.64 

0.145 

0.111 

134 

4 

4821 

4682 

2.88 

0.139 

.  0.138 

87  (all) 

Table  3.3:  Results  of  VCFEM  simulation  with  Geubelle’s  model:  (i)  the  overall  undamaged  stiffness 
Kud  for  the  2-D  section,  (ii)  the  damaged  overall  stiffness  Kd,  (iii)  the  change  in  overall  stiffness, 
(iv)  the  averaged  tensile  strain  at  the  onset  of  debonding  (v)  the  corresponding  gage  strain 
and  (vi)  the  final  debonded  angle  6^. 


fiber  # 

LAP 

NND  (mm) 

£  =  0.96 

£  =  1.04 

£  =  1.12 

Di 

Pe 

Di 

Pe 

Di 

Pe 

1 

34.26% 

0.819 

D 

D 

D 

D 

D 

D 

2 

25.75% 

2.388 

B 

B  , 

D 

D 

D 

D 

3  1 

25.76% 

0.916 

B 

B 

B 

D 

D 

D 

4 

25.12% 

1.336 

B 

B  i  B 

B 

D 

D 

5 

21.30%^ 

2.390 

B 

B  i  B 

B 

B 

D 

6 

11.01% 

6.453 

B 

B  ;  B 

B 

B 

B 

Table  3.4:  The  relation  between  quantitative  descriptors  of  the  microstructure  and  the  damage 
propagation  for  displacement  (Di)  and  periodic  (Pe)  boundary  conditions.  Values  of  strain  should 
be  read  as  e  x  D  stands  for  debonded  and  B  is  for  bonded  at  a  particular  strain. 
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1.5 


a 


*^  0  00  0.02  0.04  0  06  0.08  0.10 

Axial  Strain 


Figure  3.2:  Plot  of  the  effective  stress  Eg  as  a  function  of  axial  strain  eo  for  comparison  with  results 
in  [22]. 
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I  LoMdmg  Direction 


(a)  (b) 


Figure  3.3:  (a)  A  schematic  diagram  of  the  cruciform  specimen  with  reinforced  fibers  and  applied 
loading,  (b)  a  typical  cross-section  delineating  critical  geometric  parameters. 
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(a) 


(b) 


Figure  3.4:  (a)  Faceview  of  the  debonded  cruciform  specimen  showing  dye  penetration,  (b)  the 
cross  section  indicating  debonding  angle  as  the  limits  of  the  dye  penetrated  region. 
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Figure  3.5:  Macroscopic  pre-  and  post  debonded  stress-strain  response  of  the  specimen  with  a  single 
reinforced  circular  fiber.  The  roman  numerals  (I,  II,  III  and  IV)  correspond  to  critical  property 
values  used  to  calibrate  the  cohesive  zone  model. 
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2.0 


• - •  2-D  Plane  Strain 

»  ■2-D  Plane  Strain  a^*0.8443 

♦  2-D  Plane  Strain 
^  -  A  2-D  Plane  Strain  o„, *0.8443 


’^0.0  200  400  600  800  100.0 

Angular  Position  (degrees) 

Figure  3.6:  A  plot  of  the  normal  and  tangential  stresses  at  a  perfect  interface  by  ANSYS,  for 
obtaining  the  conversion  factor  between  3-D  and  2-D  analysis. 
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• — •  Circular  fiber.  Geubeile 
•  ■  Circular  fiber,  Needleman  Polynomial 

♦  Circular  fiber.  Needleman  Exponential 
'  ^  -4  Elliptical  fiber.  a^=0“ 

Elliptical  fiber  a^=90'’ 
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°  0.05  *  oTo  ~  ,*  9^ 

Averaged  Strain  (%) 


0.20 


Figure  3.7:  Plot  of  the  debonding  angle  6^  as  a.  function  of  average  strain  with  different  cohesive 
zone  models 


33 


^'^O.o  500  100.0  150.0  ^0.0 

Angular  position  (degrees) 

Figure  3.8:  Distribution  of  the  normal  stress  along  the  circular  interface  for  increasing  values 
of  the  debonded  angle  dj. 


Figure  3.9:  Distribution  of  the  tangential  stress  ant  along  the  circular  interface  for  increasing  values 
of  the  debonded  angle  Od- 


20.0 


Figure  3.10;  Macroscopic  pre-  and  post  debonded  stress-strain  response  of  the  specivtien  with  a 
single  reinforced  elliptical  fiber,  oriented  along  and  perpendicular  to  the  loading  axis'. 
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Figure  3.11:  Distribution  of  the  normal  stress  C7nn  along  the  elliptical  interface  for  increasing  values 
of  the  debonded  angle  9^,  {oimaj  =  O^*)- 
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6.0 


Figure  3.12:  Distribution  of  the  tangential  stress  ant  along  the  elliptical  interface  for  increasing 
values  of  the  debonded  angle  0^,  [omaj  =  O®)- 


Before  Debonding 


• . «  e=78' 

d.  -  -  o  A=00® 


Angular  position  (degrees) 

Figure  3.13:  Distribution  of  the  normal  stress  cr„„  along  the  elliptical  interface  for  increasing  values 
of  the  debonded  angle  9^,  {oimaj  ~  90°). 


39 


•3 - 0  Before  Debonding 


10.0 


'^^’^0.0  500  100.0  150  0  ^0-0 

Angular  position  (degrees) 

Figure  3.14:  Distribution  of  the  tangential  stress  ant  along  the  elliptical  interface  for  increasing 
values  of  the  debonded  angle  9d,  {ocmaj  =  90'’)- 


Gage  Strain  (%) 

Figure  3.15:  Macroscopic  pre-  and  post  debonded  stress-strain  response  of  the  specimen  with  five 
circular  fibers. 
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0  (degrees) 


Figure  3.16: 
specimen. 
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Plot  of  the  debonding  angle  0  as  a  function  of  strain,  for  each  fiber  in  the  multi- fiber 
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Figure  3.17:  Tessellated  mesh  of  Voronoi  elements  for  microstructures  containing  100  fibers  with 
area  fraction=20%,  (a)  Uniform  distribution  (b)  Random  distribution. 


(b) 

Figure  3.18:  (a)  Macroscopic  pre-  and  post  debonded  stress-strain  response  of  the  microstructures 
containing  100  circular  fibers  and  (b)  microscopic  plot  of  the  number  of  debonded  fibers  as  a  func¬ 
tion  of  strain.  The  6  plots  in  each  figure  correspond  to:  (1)  vf=20%,  uniform  dispersion-uniform 
size-periodic  boundary,  (2)  vf=20%,  uniform  dispersion-uniform  size-displacement  boundary,  (3) 
vf=20%,  random  dispersion-uniform  size-displacfjnent  boundary,  (4)  vf=20%,  random  dispersion- 
uniform  size-periodic  boundary,  (5)  vf=10%,  random  dispersion-uniform  size-displacement  bound¬ 
ary  and  (6)  vf=10%,  random  dispersion-random  size-displacement  boundary. 


Max.- 


Min.  - 


-  1.495E-01 
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8.970E-02 
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-  O.OOOE+00 


(a) 


Max.  - 
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5.928E*02 


2.964E-02 


0,0O0E+O0 


(b) 

Figure  3.19:  Contour  plots  of  the  microscopic  axial  stros.s  a for  the  unformly  dispersed  microstruc¬ 
ture  (vf=20%)  at  €  =  1.34  X  10"‘‘;  (a)  without  periodicity  and  (b)  with  periodicity. 
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Figure  3.20:  Contour  plots  of  the  micros 
crostructure  (vf=20%)  at  e  =  1,12  x  10“"^; 
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Figure  3.21:  Distribution  of  the  radial  stress  Orr  along  the  interface  for  different  fibers  (1,  2,  3 
and  4  in  figure  3.17a)  of  the  uniform  raicrostructure  with  displacement  boundary  conditions  for  (a) 
e  =  1.31  X  10"''  and  (b)  e  =  1.33  x  10”'* 
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t).!  0.15  02  0.25  0.3  0.35  0.4 

Local  area  fraction 

(a) 


Nearest  neighbor  S-S  distance 

(b) 

Figure  3.22:  Histograms  showing  the  number  of  debonded  and  bonded  fibers  as  functions  of  (a) 
local  area  fraction  (LAF)  (b)  nearest  neighbor  surface  to  surface  distance  (NND),  with  periodicity 
boundary  conditions.  1,  2,  3  in  the  legend  correspond  to  increasing  strain  levels  e  =  0.96  x  10“'*, 
e  =  1.04  X  lO""*  and  e  =  1.12  x  lO"**,  while  4  in  the  legend  corresponds  to  undamaged  interfaces. 
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Nearest  neighbor  S-S  distance 

(b) 

Figure  3.23;  Histograms  showing  the  number  of  deboiuled  and  bonded  fibers  as  functions  of  (a)  local 
area  fraction  (LAF)  (b)  nearest  neighbor  surface  to  surface  dismnce  (NND),  with  displacement 
boundaury  conditions.  1,  2,  3  in  the  legend  correspond  to  increasing  strain  levels  e  =  0.96  x  10“'*, 
€  =  1.04  X  lO"'*  and  €  =  1.12  x  10“*,  while  4  in  the  legend  corresponds  to  undamaged  interfaces. 
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Chapter  4 


A  Multi-level  Computational  Model 
for  Multi-scale  Damage  Analysis  in 
Composite  and  Porous  Materials 


4.1  Introduction 

Heterogeneous  structures  with  second  phase  inclusions  or  voids  in  the  microstructure  are  conven¬ 
tionally  analyzed  with  macroscopic  properties  obtained  from  homogenization  of  response  at  smaller 
(meso-.  micro-)  length  scales.  The  mathematical  homogenization  theory,  which  uses  asymptotic 
expansions  of  displacement,  strain  and  stress  fields  about  macroscopic  values,  has  been  used  as  a 
tool  for  analyzing  multiple  scale  responses  in  Benssousan  et.  al.  (1978),  Sanchez-Palencia  (1980), 
Parton  and  Kudryavtsev  (1993),  Bakhvalov  and  Panasenko  (1984).  The  method  is  based  on  as¬ 
sumptions  of  spatial  periodicity  of  microscopic  representative  volume  elements  (RVE)  and  local 
uniformity  of  macroscopic  fields  within  each  RVE.  It  decomposes  the  multiscale  boundary  value 
problem  into  a  decoupled  set  of  micro-scale  RVE  problem  and  a  macro-scale  problem.  Concurrent 
finite  element  analyses  are  executed  at  the  each  scale  for  information  transfer  between  the  scales. 
.Multiple  scale  analysis  of  linear  elastic  reinforced  composites  by  this  method  have  been  conducted 
by  Fish  and  Belsky  (1995),  Fish' and  Wagiman  (1993).  Guedes  and  Kikuchi  (1991)  and  Hollister 
and  Kikuchi  (1992).  For  nonlinear  materials,  the  homogenization  methods  have  been  extended  by 
Suquet  (1985),  Fish  et.al.  (1997),  Guedes  (1990)  and  Cheng  (1992).  The  method  has  also  been 
implemented  to  simulate  damage  by  fiber-matrix  debonding  in  linear  elastostatics  by  Lene  (1986) 
and  fiber  rupture  using  a  phenomenological  damage  model  by  Devries  et.  al.  (1989). 

Despite  its  advantages,  asymptotic  homogenization  Inus  suffered  shortcomings  arising  from  effi¬ 
ciency  and  accuracy  considerations.  Enormous  computational  efforts  can  result  wth  this  method 
due  to  the  fact  that  at  each  integration  point  in  the  macroscopic  model,  boundary  value  problems 
of  the  microstructural  RVE  should  be  solved  twice.  To  ironomize  computations,  many  studies 
have  assumed  simple  unit  cells  models  of  the  microscopic  RVE.  Such  idealizations  may  however  be 
unrealistic  for  deformation  and  failure  analysis  of  many  materials.  The  homogenization  method 
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has  another  major  limitation  stemming  from  its  basic  assumptions,  viz.  (a)  uniformity  of  the 
macroscopic  fields  within  each  RVE  and  (b)  spatial  periodicity  of  the  RVE.  The  uniformity  as¬ 
sumption  is  not  appropriate  in  critical  regions  of  high  gradients,  where  the  macroscopic  helds  can 
yg^y  considerably.  Free  edges,  interfaces,  macrocracks,  neighborhood  of  material  discontinuities 
and  most  importantly  in  the  regions  of  evolving  microscopic  damage  and  instability  are  potential 
sites  of  nonuniformity.  Furthermore,  statistical  periodicity  implies  that  the  RVE  may  be  repeated 
to  represent  the  entire  neighborhood  of  a  macroscopic  point.  For  non-uniform  microstructures  suf¬ 
ficiently  large  portions  should  be  considered  as  RVE  for  homogenization  analysis.  Unit  cell  models 
are  severely  limited  in  this  respect.  Even  higher  order  theories  of  homogenization  may  be  com¬ 
putationally  unviable.  A  few  effective  global-local  techniques  based  on  hierarchical  decomposition 
and  superposition  of  field  variables  have  been  proposed  by  Belytschko  et.  al.  (1994),  Robbins 
and  Reddy  (1996)  and  Hughes  (1995).  Pagano  and  Rybicki  (1974)  had  discussed  the  breakdown 
of  effective  modulus  theory  for  composite  laminates  with  free  edges  and  the  need  for  global-local 
techniques.  Fish  and  Belsky  (1995)  and  Fish  et.  al.  (1994)  have  used  global-local  techniques  with 
multigrid  methods  to  extend  the  multiple  scale  modeling  to  non-periodic  materials.  Zohdi  et.  al. 
(1996),  Oden  and  Zohdi  (1997)  and  Oden  et.  al.  (1999)  have  developed  a  homogenized  Dirichlet 
projection  method  (HDPM),  which  resolves  the  microstructural  effects  at  different  scales  on  the 
macroscopic  response  of  heterogeneous  structures. 

Adaptivity  in  the  computational  modules  for  multiple  scale  problems  entails  minimizing  two 
types  of  errors,  viz.  the  discretization  error  and  the  modeling  error  due  to  homogenized  material 
properties  as  discussed  in  Oden  et.  al.  (1997,1999).  In  the  present  work,  an  adaptive  multilevel 
method  is  proposed  in  to  address  the  latter  type  of  modeling  error.  It  is  aimed  at  improving 
the  accuracy  of  analyses  of  elastic-plastic  composite  and  porous  structures  with  microstructural 
damage.  The  model  uses  computational  hierarchy  to  concurrently  predict  the  evolution  of  vari¬ 
ables  at  structural  aind  microstructural  scales,  as  well  as  to  track  the  incidence  and  propagation  of 
microstructural  damage.  Analysis  of  microstructural  response  with  arbitrary  distributions,  shapes 
and  sizes  of  heterogeneities  is  conveniently  done  by  the  Voronoi  Cell  finite  element  model  (VCFEM) 
(seee.g.  Moorthy  and  Ghosh  (1996,1998),  Ghosh  et.  al.  (1997)).  A  high  level  of  computational 
efficiency  with  sufficient  accuracy  and  resolution  has  been,  achieved  for  elastic  and  elastic-plastic 
materials  by  this  method.  Progressive  damage  by  particle  cracking  has  been  done  in  Moorthy 
and  Ghosh  (1998).  Recently,  an  adaptive  VCFEM  has  been  successfully  proposed  by  Moorthy  and 
Ghosh  (1999),  where  element  adaptation  is  executed  in  two  consecutive  stages  based  on  a-posteriori 
evaluation  of  error  measures.  In  the  first  stage,  displacement  function  adaptations  are  carried  out 
by  a  /i-refinement  and  p-enrichment  strategy,  which  is  followed  by  an  enrichment  of  stress  func¬ 
tions  to  reduce  the  error  in  kinematic  relations. 

For  periodic  representative  volume  elements  (RVE)  of  elastic  and  elastic-plastic  materials,  the 
microstructural  VCFEM  has  been  coupled  with  structural  analysis  codes  by  qsing  asymptotic 
homogenization  in  Ghosh,  Lee  and  Moorthy  (1995,1996).  This  method  fails  for  problems  where  as¬ 
sumptions  of  m2w:roscopic  uniformity  and  statistical  periodicity  are  questionable.  Consequently,  it 
becomes  necessary  to  implement  a  combination  of  homogenization  and  global-local  methods,  which 
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is  a  nontrivial  undertaking  due  to  lack  of  apriori  knowledge  of  regions  that  require  differential 
resolution.  The  multi-level  methodology  developed  in  this  work  addresses  this  issue  by  differenti¬ 
ating  between  non-critical  and  critical  regions  and  ranging  from  macroscopic  computations  using 
continuum  constitutive  relations  to  zooming  in  at  ‘hotspots’  for  pure  microscopic  simulations.  The 
zoom-in  is  accomplished  by  a  mesh-enrichment  technique  or  /i-adaptation,  where  macroscopic  el¬ 
ements  are  successively  dissected  in  regions  of  steep  solution  gradients.  Also  the  adaptations  are 
based  on  local  estimates  of  ‘error’  or  solution  gradients.  The  work  introduces  three  levels  of  the 
computational  domain,  for  (a)  fully  macroscopic  analysis  with  homogenized  material  parameters, 
(b)  macro-microscopic  analysis  for  periodic  RVE’s  and  (c)  fully  microscopic  analysis.'  Additionally 
a  new  way  of  developing  a  piecewise  continuous  elastic-plastic  constitutive  model  is  presented  in 
this  work.  This  model  accounts  for  the  details  of  microstructural  morphology  and  variables.  A 
number  of  numerical  examples  are  solved  with  vjurious  microscopic  architectures  to  support  the 
development  of  this  multi-level  model. 

4.2  Two  Way  Coupling  for  Multiple  Scale  Modeling 

Multiple  scale  modeling  of  heterogeneous  materials  is  necessary  to  concurrently  account  for  spatial 
variability  at  the  macro-  and  micro-scales.  An  effective  model  of  this  class  requires  two-way  coupling 
for  efficient  computing,  as  well  as  for  accurate  representation  of  the  necessary  variables  at  different 
scales.  The  first  is  a  'bottom  up’  coupling  for  determination  of  equivalent  homogeneous  behavior  at  a 
macroscopic  point  x,  as  a  function  of  the  microstructural  geometry  and  behavior  of  the  constitutive 
phases,  but  independent  of  applied  loads  to  the  structure.  In  the  homogenization  procedure,  an 
isolated  representative  volume  element  (RVE)  Y{x)  C  3?^  is  identified  at  microstructural  scale 
of  heterogeneities  (figure  ??b).  The  scale  y  of  the  RVE  domain  y(x)  is  large  with  respect  to 
the  characteristic  length  1  of  microscopic  heterogeneities,  but  is  significantly  small  compared  to 
the  macroscopic  length  scale  L  of  the  structure  and  applied  loads.  Homogenized  variables  at  the 
macroscopic  scale  are  obtained  by  volume  averaging  of  variables  in  the  RVE,  following  the  definition 

<  /  >K(x)=  fiy)dy  where  \Y\  =  j^dY  (4.1) 

The  condition  for  macroscopic  homogeneity,  according  to  the  Hill-Mandel  hypothesis  (see  e.g.  Hill 
(1965)),  assumes  equivalence  of  strain  energy  for  the  actual  and  equivalent  homogenized  media. 
Thus  for  a  statically  admissible  stress  field  a[y)  and  kinem.atically  admissible  strain  field  e(y), 

<  <7  :  €  >y=<  <7  >y:<  e  >v  Vy  €  T  (4.2) 

The  microscopic  stress  <r{y)  and  strain  e(y)  fields  satisfying  the  homogeneity  condition  (4.2)  may 
be  obtained  by  solving  boundary  value  problems  for  the  RVE  Y  with  prescribed  homogeneous 
stress /strain  or  periodicity  boundary  conditions,  stated  as:  e. 

T**  =<  (7  >y  •n(y)  =  <7  ■  n(y)  on  dY  :  Uniform  Traction  (a) 

u**  =<  €  >y  y  on  dY  :  Uniform  Strain  (b) 


51 


£  >y  -y  +  u(y)  =<  c  >/  y  +  u(y  +  kY)  on  dY  :  Y-Periodicity  (c)  (4.3) 

where  k  is  an  3x3  array  of  integers  and  Y  is  the  period  of  V-periodic  displacement  functions  u, 
interpreted  as  local  perturbations  to  macroscopic  strain  based  displacement  fields.  The  macroscopic 
constitutive  equations  are  obtained  by  solving  a  boundary  value  problem  of  the  RVE  Y  with  any 
one  of  the  three  sets  of  boundary  conditions  in  equation  (4.3),  followed  by  the  averaging  process  in 
equation  (4.1).  For  linear  elastic  constituent  phases  in  Y,  the  relation  between  the  strain  energy 
functions  has  been  established  in  Suquet  (1987)  as: 

<  e  >;  :<  e  >  <  <  e  >:  :<  e  >  <  <  €  >:  >  VE  C  3?®  (4.4) 

where  E{)..  E\r-  Eitr  are  respectively  the  homogenized  stiffness  tensors  evaluated  with  uniform  trac¬ 
tion,  periodicity  and  uniform  strain  boundary  conditions,  and  <  e  >  is  the  macroscopic  (applied 
or  averaged)  strain  field.  The  difference  in  stiffnesses  with  the  kinematic  and  kinetic  boundary 
conditions,  reduce  with  decreasing  size  of  Y.  It  is  generally  concluded  by  Hollister  and  Kikuchi 
(1992)  that  for  the  same  RVE  size,  the  periodicity  boundary  conditions  are  expected  to  yield  more 
accurate  statistically  homogenized  constitutive  parameters  and  macroscopic  properties. 


The  other  coupling  is  the  ‘top  down’  where  the  evolution  of  variables  are  evaluated  in  the 
microstructure  from  known  macroscopic  variables,  by  a  process  termed  as  localization.  In  those 
regions,  where  the  microstructure  may  admit  a  RVE  Y ,  the  microscopic  variables  can  be  evaluated 
by  solving  a  boundary  value  problem  with  imposed  macroscopic  strains  and  the  local  periodicity 
condition  in  equation  4.3c.  In  other  regions,  where  assumptions  of  local  periodicity  of  the  RVE 
may  be  unrealistic,  the  localization  process  will  entail  direct  interfacing  of  the  microstructural  and 
macroscopic  regions. 

4.3  A  Multi  Level  Model  for  Coupling  Different  Scales 

In  the  spirit  of  true  two  way  coupling  of  multiple  scale  problems,  the  computational  domain  in  this 
work  is  adaptively  decomposed  into  three  levels  of  hierarchy  based  on  requirements  of  resolution. 
Such  hierarchy  is  intended  to  increase  computational  efficiency  as  well  as  accuracy  in  concurrent 
prediction  of  variables  at  the  continuum  and  microstructural  scales.  As  proposed  in  Lee  et.  al. 
(1999),  the  model  uses  homogenization  of  microstructural  RVE  solutions  to  evaluate  homogenized 
properties  and  cascades  down  the  scales  at  hotspots  of  evolving  damage.  The  three  levels  of  hier¬ 
archy  with  requirements  of  increasing  resolution  (figure  4.2)  are  as  follows. 


•i.  Computational  Subdomain  Level-0; 

These  correspond  to  non-critical  macroscopic  regions  in  figure  4.2a.  where  deformation  variables 
are  relatively  uniform  and  periodicity  conditions  may  be  assumed  for  the  underlying  material  RVE. 
Scale  effects  are  negligible  in  this  region  and  local  constitutive  relations  may  be  derived  from  pos¬ 
tulates  of  the  RVE  approaching  zero  volume.  Continuum  level  anisotropic  plasticity  constitutive 
relations,  that  are  consistent  with  the  actual  microstructural  constitution,  are  developed  for  macro¬ 
scopic  modeling  of  these  regions. 
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The  level-0  macroscopic  simulations  are  accompanied  by  element  refinement  or  h- adaptation  for 
two  reasons.  The  first  is  to  identify  and  reduce  a  chosen  ‘error  measure’  in  the  macroscopic  com¬ 
putational  model.  A  second  attribute  is  that  it  enables  the  computational  model  to  zoom  in  on 
regions  of  evolving  nonuniformity  due  to  microscopic  non- homogeneity.  This  reduces  the  disparity 
in  size  of  the  macro-  and  micro-  scale  elements  by  successive  refinement  of  macroscopic  elements 
in  the  critical  regions  as  shown  in  figure  4.2a. 

•ii.  Computational  Subdomain  Level-1:  These  are  regions  that  face  imminent  microscopic 
non-homogeneity  and  resulting  macroscopic  nonuniformities  (figure  4.2a, b).  Though  the  computa¬ 
tions  are  still  macroscopic,  concurrent  monitoring  of  the  development  of  damage  and  instabilities  in 
the  RVE  is  possible  in  this  level.  For  concurrent  macro-  and  micro-  scale  analyses  the  asymptotic 
homogenization  methods,  which  is  based  on  the  existence  of  an  RVE,  is  used.  Macroscopic  element 
refinement  by  h-adaptation  continues  for  this  level. 

#iii.  Computational  Subdomain  LeveI-2;  These  critical  regions  materialize  with  the  evolution 
of  microstructural^ damage  in  the  form  of  evolved  microcracks  or  instabilities  (figure  4.2b),  leading 
to  high  macroscopic  field  gradients.  The  assumptions  of  macroscopic  uniformity  and  local  periodic¬ 
ity  are  unrealistic.  To  realize  scale  effects,  it  is  required  that  the  level-1  macro-micro  computational 
model  switch  to  a  completely  microscopic  model,  encompassing  large  portions  of  the  microstruc¬ 
ture.  A  detailed  flow  chart  of  the  adaptive  hierarchical  process  is  depicted  in  figure  4.3. 

The  substructured  computational  domain  is  delineated  as  an  elastic-plastic  body  of  material 
domain  Clm,\T  that  consists  of  regions  fo*"  which  the  RVE  is  repeated  periodically,  and 

also  of  regions  ^MAT(np)  where  the  periodicity  assumptions  do  not  hold,  i.e. 

f^A/Ar(x,  u)  =  ^MAT{p)(^^  u)  U  ^S{AT{np){^^  u)  (4.5) 

The  macroscopic  regions  of  periodicity  are  further  constituted  of  repeating  a  large  number  of  RVE’s 

''''  =  uS*?"''’''”’  where  =  uS.n'  (4.6) 

Here  VCp)  corresponds  to  the  number  of  different  RVE’s  in  the  periodic  regions,  and  for  all 
practical  purposes,  oo  corresponds  to  a  sufficiently  large  number.  The  non-periodic  region  ^MAT{np) 
is  defined  as  the  set  of  all  microstructural  regions  for  which  the  N{np)  RVE’s  are  not  repeated,  i.e. 

:  li  n  r,  W  0  V  t  (  (4.7) 

The  level-0  and  level-1  of  the  computational  domain  <>  correspond  to  the  periodic  regions,  while 
the  level-2  belongs  to  the  non-periodic  regions  as 

fljoUfi/l  C  ^MAT{p)  •  ^  ^.UAT{np) 


53 


4.4  Homogenization  with  Voronoi  Cell  FEM 


4.4.1  Asymptotic  homogenization  with  microstructural  periodicity 

Consider  a  heterogeneous  structure  occupying  a  region  ^structure  (figure  ??a),  for  which  the  hetero¬ 
geneous  microstructure  constitutes  of  spatially  repeated  RVE’s  Y (x)  about  a  macroscopic  point  x 
as  shown  in  figure  ??b.  The  RVE  is  discretized  into  a  mesh  of  Voronoi  cells,  which  naturally  evolve 
from  the  microstructure  by  Dirichlet  tessellation.  In  the  Voronoi  cell  FEM  (VCFEM).  each  Voronoi 
cell  represents  a  basic  structure'  element  (BSE)  which  is  the  neighborhood  of  a  heterogeneity  in 
the  microstructure.  The  dimensions  of  the  Y{x)  are  typically  very  small  in  comparison  with  the 
structural  dimensions  L.  i.e.  j  'S  a  very  small  positive  number  e.  Due  to  variation  of  evolutionary 
variables  in  a  small  neighborhood  e  of  the  macroscopic  point  x,  all  variables  are  assumed  to  exhibit 
dependence  on  both  length  scales  i.e.  =  4>(x,  f ),  where  y  =  7.  The  superscript  e  denotes 
association  of  the  function  with  the  two  length  scales  and  hence  fl'  corresponds  to  a  connected 
structural  and  microstructural  domain.  The  assumption  of  periodic  repetition  of  the  microstruc¬ 
ture  about  X  makes  the  dependence  of  the  function  on  y(=  7).  V-periodic  (see  e.g.  Bakhvalov  and 
Panasenko  (1984),  Guedes  and  Kikuchi  (1991),  Hollister  and  Kikuchi  (1992),  Devries  (1989)).  For 
small  deformation  elasto-plasticity,  the  rate  forms  of  the  equilibrium,  kinematic  and  constitutive 
relations  are  given  as: 

if  dill  du\ 

EquilibriuTTi  :  ^ij,j  ~  /•>  Kincmcitics  .  —  I  Qj.e  dx^ 

Constitutive  :  =  E^jf-ieh  in  (4.9) 

where  (x,y),  e'^  (x,y)  and  u'(x,y)  are  V-periodic  rates  of  stress,  strain  and  displacement  fields 
respectively.  Furthermore  the  periodic  boundary  conditions  may  be  specified  as 


u‘(x,y)  =  U'(x,y  -t-  kY)  on  RVE  boundary  dY  (4.10) 

and  olj  is  continuous  across  dY 


The  V-periodic  displacement  rate  field  is  approximated  by  an  asymptotic  expansion  about  x  with 
respect  to  the  parameter  e: 


u'(x)  =  txf(x,y)  -H€u,'(x,y)  -H€-'uf(x.y)  • 


2, -,2/ 


X 

y  =  7 


(4.11) 


Noting  that  the  spatial  x‘  derivative  of  any  function  depends  on  the  two  length  scales  and  is  given 

I  ^ 

f  dyt 


as: 


_d 

dx 


/  ,  x,\  (>4>  I  0^ 

;(t(x,y  =  -)]  =  - 


the  stress  rate  tensor  can  be  expressed  as 


(4.12) 


where 


^0  -  F‘  ^ 
<7ij  -  Ex]kl  Qy^ 


^ij  —  ^ijkl 


'dul  ^ 

dxi  dyi 


—  ^ijkl 


'duj  dill' 
dxi  dyi 


(4.14) 


From  equations  4.9  and  4.13,  and  using  the  periodicity  condition  on  the  RVE  boundary  f  dTr  = 
0,  it  can  be  proved  (see  e.g.  Ghosh  et.  al.  (1995,1996))  that 


0,  u^  =  u^(x}  and 


dyj  dyj  dyi  dyi 


(4.15) 


By  neglecting  the  terms  associated  with  e  or  higher  in  equation  (4.13),  the  constitutive  relation  in 
the  y-domain  is  expressed  as 


'V 


=  cr, 


1  — 


Eijkl^kl  —  ^ijkii 


dyi’ 


+ 


(4.16) 


Here  e%i  is  the  microstructural  strain  rate  tensor,  for  which  e*,/  =  ^(^  +  5^)  is  an  averaged 

macroscopic  part  and  ^ +  5^)  is  denoted  as  a  fluctuating  strain  rate  tensor  (see  e.g. 

Suquet  (1985)).  Due  to  linearity  of  the  rate  problem,  &lj,  ii}  and  the  microscopic  equilibrium 
condition  can  be  expressed  as 


_  :^klj  \^^k 

<^xj  -  (y)  ’ 


u}  =  X\ 


dxi 


ddfj(y) 

dyj 


=  0 


(4.17) 


In  equation  (4.17),  dfj  is  a  F-antiperiodic  function  and  xf  ‘S  a  F-periodic  function  representing 
characteristic  modes  of  the  deformation  in  the  RVE.  Substituting  equation  (4.17)  in  the  constitutive 
relations  of  equation  (4.9)  yields  the  microscopic  constitutive  relations  as: 


d^j(y)  =  ^ijpml^kp'^lm  + 


dx‘: 


kl 


dyn 


(4.18) 


where  (5,j  is  Kronecker  delta.  The  mean  of  equation  (4  18)  yields  the  homogenized  elastic-plastic 
tangent  modulus  for  use  in  the  macroscopic  analysis,  (n  the  form 


^ijkl 


=  < 


A*' 


>Y 


-u 


df,dY  = 


1 

1Y| 


dXp, 


^ijpmi^kp^lm  +  )dY 


(4.19) 


The  macroscopic  stress  and  strain  relation  can  thus  be  stated  as 

a^mn  n,-0 

Eij(x)  =<  E^ji^iiSkm^ln  H - >y—  ^ijmn^rnni'^)  (4.20) 

where  the  homogenized  variables  are  S(x)  =<  cr'(x,y)  >y  and  e(x)  =<  e*fx,y)  >y-  The 
incremental  small  deformation  analysis  for  elastic-plastic  materials  can  be  conducted  with  the 
homogenized  modulus  at  the  macroscopic  level  and  by  using  the  Voronoi  cell  finite  element  model 
(VCFEM)  for  solving  the  microscopic  problem. 
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4.4.2  The  Voronoi  Cell  FEM  for  Microstructural  Analysis 

The  Voronoi  cell  finite  element  model  (VCFEM)  has  been  successfully  developed  for  composite  and 
porous  materials  in  Moorthy  and  Ghosh  (1996,1998,1999).  Arbitrary  dispersion  patterns,  shapes 
and  sizes  of  heterogeneities  are  readily  modeled  by  VCFEM.  The  computational  model  naturally 
evolves  by  Dirichlet  tessellation  of  the  microstructure  as  shown  in  figure  ??b.  Each  Voronoi  cell 
with  the  embedded  inclusion  or  void  is  treated  as  an  element  in  this  formulation.  In  Moorthy 
and  Ghosh  (1998)  and  Ghosh  and  Moorthy  (1998),  the  VCFEM  formulation  has  been  extended 
to  include  damage  evolution  in  the  form  of  particle  cracking,  where  the  crack  is  realized  as  an 
elliptical  void.  Each  Voronoi  cell  element  is  amenable  to  change  in  topology  frorr  iwo  constituent 
phases  (matrix  and  inclusions)  in  undamaged  cells,  to  three  phases  (matrix,  inclusion  and  crack) 
in  damaged  cells.  Complete  particle  cracking  or  splitting  is  assumed  to  occur  at  the  very  onset  of 
damage. 


The  VCFE  formulation  constructs  a  hybrid  element  by  combining  the  aspects  of  finite  element 
methods  with  important  micromechanics  considerations.  Use  of  a  hybrid  stress  based  formulation 
results  in  a  high  level  of  accuracy  with  a  significantly  reduced  degree  of  freedom,  compared  to 
displacement  based  FEM  models.  Consider  a  typical  representative  volume  element  Y  consisting  of 
.V  undamaged  and/or  damaged  particles,  that  are  contained  in  each  of  the  N  Voronoi  cell  elements, 
as  shown  in  figure  ??(b).  In  VCFEM,  the  RVE  Y  is  comprised  of  the  Voronoi  cell  elements  Ye,  i.e. 
Y  =  ^e=\^e-  The  assumed  stress  hybrid  formulation  in  VCFEM  requires  independent  assumptions 
of  an  equilibriated  stress  field  (ff)  in  each  of  the  matrix  and  inclusion  phases  of  each  element 
=  y"i  y^c^  compatible  displacement  fields  u  on  the  element  boundary  dYe,  u'  on  the 
matrix-inclusion  interface  dY^  and  on  the  crack  boundary  dYcr-  In  3'^  incremental  formulation  for 
elasto-plasticity,  the  incremental  variational  formulation  introduces  an  element  energy  functional. 


nf  (Act,  Au)  =  -  f  AB{(r,  Act)  dY  -  f  e  :  Aa  dY 

JYc  •'Ye 

-h  /  (cr  -f  A<t)  ‘  *  (u  +  Au)  ddY  (Inter-element  Traction  Reciprocity) 

JdYe 

-  [  (t  -f  At)  •  (u  +  Au)  dV  (Boundary  Traction) 

JVtm 

-  f  ((T^  -h  ‘  ♦  (u'  +  Au')  dY  (Matrix-Inclusion  Interface  Traction) 

JdYc 

-  /  (cr^  -h  Aa^)  •  •  (u"  +  Au")  dY  (Crack  Boundary  Traction)  (4.21) 

JdYrr 


where  AB  is  the  increment  of  complimentary  energy  density.  Variables  (<t,u)  correspond  to  values 
at  the  beginning  of  an  increment,  while  variables  (Act,  Au)  are  the  corresponding  increments  in 
a  load  increment  or  step.  Outward  normals  on  dYg,  dYc  and  dYcr  denoted^ by  n  ,  n  and 
respectively.  Superscripts  ttt,  c  and  cr  are  associated  with  the  matrix,  inclusion  and  crack 
phases  respectively  in  each  Voronoi  cell  element.  The  total  energy  for  the  entire  RVE  of  N  Voronoi 
cells  is  obtained  as  11^  =  nf .  Setting  the  first  variation  of  Ilf  in  equation  (4.21)  with 

respect  to  stress  increments  Act  to  zero  yields  the  element  compatibility  as  the  Euler  equation, 


56 


while  setting  the  first  variations  of  with  respect  to  the  independent  boundary  displacements 
Au,  Au'  and  Au"  to  zero,  yield  the  inter-element  boundary  traction  reciprocity,  matrix-inclusion 
interface  traction  reciprocity  and  zero  traction  crack  boundary  condition  respectively.  Independent 
assumptions  on  stress  increments  Atr  are  made  in  the  matrix  and  inclusion  phases  in  each  element, 
thus  allowing  stress  discontinuities  across  the  interface.  In  this  process  special  forms  of  the  Airy’s 
stress  function  $(x,y)  to  enhance  computational  efficiency,  has  been  developed  in  Moorthy  and 
Ghosh  (1996.1998)  for  equilibriated  stress  fields.  The  functions  facilitate  stress  concentration  near 
the  interface  and  crack  boundary,  accounting  for  the  shape  of  the  inclusion  and  crack  and  also 
help  satisfy  traction  reciprocity  at  the  interfaces  dYc  s^nd  dYcj-.  Furthermore,  they  decay  at  large 
distances  from  the  interfaces.  Compatible  displacement  increments  are  generated  on  each  of  the 
boundaries/ interfaces  5^,  dYc  and  dYcr  by  interpolating  nodal  displacements  using  polynomial 
shape  functions.  The  stress  and  displacement  interpolations  may  be  expressed  as: 

{Act'"}  =  [P'"(i,y)]{Aj3'"}  (in  matrix)  and  (Act"}  =  [P^(x,y)]{A/9‘'}  (in  inclusion) 

(Au)  =  [L®]{Aq}  (on  element  boundary)  .  (Au'}  =  [L‘=]{Aq'}  (  on  interface),  and 
(Au"}  =  [L‘^’‘]{Aq"}  (  on  crack  face)  (4.22) 

where  (Aq).  (Aq'}  and  (Aq"}  are  the  nodal  displacement  increment  vectors,  and  [L^],  [L'’]  and 
are  the  corresponding  interpolation  matrices. 

4.4.3  Coupling  Asymptotic  Homogenization  with  VCFEM 

In  the  incremental  formulation,  the  equilibriated  microscopic  stress  increment  corresponds  to 
Acr^(=  Act')  in  equation  (4.17)  and  the  microstructural  strain  increments  are  designated  as  Ae'  in 
equation  (4.16).  Similarly,  the  increments  in  microscopic  displacements  on  the  cell  boundaries  dYg 
are  identified  with  Au^  in  equation  (4.17)  and  those  on  the  interface  and  crack  surface  are  denoted 
by  Au*’  and  Au*”  respectively.  In  the  absence  of  explicit  traction  boundaries  due  to  periodicity 
conditions  on  the  boundary,  the  incremental  energy  functional  for  each  Voronoi  cell  element  in 
equation  (4.21)  is  modified  for  the  asymptotic  homogenization  process  as: 


nf  =  -  €<  :  A<r<  dY  + 

f  (cT'^AcT)'  n'-(u*+Au*)day  -  [  (cr''"  +  Act'’"  -  cr'"  -  Act'")  •  n"  •  (u*' -b  Au*')  97 

JdY, 

-  [  (ct'"  +  Act'")  •  n"'"  •  (u*"  -b  Au*")  dY  -r  [  (e  +  Ae)  Acr'c^F  (4.23) 

JdYcr 

where  is  an  instantaneous  elastic-plastic  compliance  tensor.  The  last  term  in  equation  (4.23) 
incorporates  the  effect  of  macroscopic  strains  e  in  the  microstructure.  The  stationary  condition 
of  nf  with  respect  to  stress  increment  Acr'j  yields  as  Euler's  equations,  the  inctemental  form  of 
kinematic  relations 


e,j  -b  Ae,j  e,j 


^(^t  "b  jjjg  element  matrix  and  inclusion  phases  (4.24) 
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Stationaxity  of  the  total  energy  functional  Uf  with  respect  to  displacement  increments 

Au- ,  Au-'  and  AuJ”,  result  in  the  inter-element,  interface  and  crack  boundary  traction  reciprocity 
conditions  respectively. 

(4  +  A4)-nf  =  -(4+Aa'j)-nf  on  dYe 
(<t‘; -h  A4")  •  +  Aa'f)  •  on  dYc 

{(T^  +  Aaf/)  •  nf  =  0  on  dYcr  (4.25) 

The  microscopic  VCFEM  module  is  executed  for  two  purposes  in  each  increment  of  the  macro¬ 
scopic  module.  The  first  is  to  evaluate  the  microscopic  stress  increments  A<t^  from  given  values 
of  the  macroscopic  strain  e  at  the  beginning  of  the  step,  and  its  increment  Ae.  The  second  is  to 
calculate  the  instantaneous  homogenized  tangent  modulus  at  the  end  of  the  increment  in  the 
macroscopic  module.  The  details  of  the  calculation  of  microscopic  stress  and  homogenized  tangent 
modulus  calculation  procedure  are  given  in  Ghosh  et.  al.  (1995,1996).  The  Voronoi  cell  finite  ele¬ 
ment  module  is  incorporated  in  a  macroscopic  analysis  module  with  the  interface  being  created  by 
the  homogenization  procedure.  The  macroscopic  analysis  is  performed  using  a  displacement  based 
finite  element  code  with  plane  strain  QUAD4  elements.  Numerical  integration  in  these  elements 
uses  one-point  reduced  integration  and  hourglass  control  (see  Koh  and  Kikuchi  (1987)).  Material 
constitutive  relations  at  each  integration  point  of  elements  are  obtained  from  homogenization  of 
microscopic  VCFEM  results.  Microscopic  stresses  af^are  averaged  to  yield  macroscopic  stresses 
E.j.  The  microscopic  VCFEM  is  also  invoked  to  evaluate  the  homogenized  elastic-plastic  tangent 
modulus  by  applying  unit  components  of  macroscopic  strain. 

4.5  Computational  Subdomain  Level-0  in  the  Hierarchical  Model 

Level-0  corresponds  to  macroscopic  regions  (fi/o  ^  ^xfAT{p))  figure  4.2a,  where  deformation  vari¬ 
ables  like  stresses  and  strains  are  relatively  uniform  in  their  macroscopic  behavior.  Scale  effects  are 
negligible  and  local  constitutive  relations  may  be  derived  from  postulates  of  zero  RVE  volume  in 
the  limit  and  periodicity.  It  is  assumed  that  macroscopic  analyses  with  homogenized  constitutive 
relatio.iS  are  sufficient  for  these  regions.  Anisotropic  constitutive  relations  with  varying  param¬ 
eters  are  developed  for  continuum  analysis  of  heterogeneous  microstructures  with  elastic-plastic 
constituent  phases.  To  account  for  details  of  microstructural  morphology,  the  constitutive  model 
is  based  on  two  scale  analysis  using  the  asymptotic  homogenization  method  and  microstructural 
analysis  by  VCFEM.  A  continuum  constitutive  model  can  greatly  enhance  computational  efficiency 
over  twoscale  analysis. 


4.5.1  An  Elastic-Plastic  Constitutive  Model 

V^arious  continuum  constitutive  models  have  been  proposed,  based  on  unit  cell  analyses  of  com¬ 
posite  and  porous  microstructures.  One  parameter  plastic  potential  functions  with  assumptions  of 
anisotropy  have  been  introduced  in  Sun  and  Chen  (1991)  and  Xie  and  Adams  (1995)  for  composite 
materials,  where  the  parameter  is  determined  by  least  squares  fitting  of  unit  cell  characteristic 


58 


responses.  Bao  et.  al.  (1991)  have  used  the  same  hardening  exponent  for  the  composite  as  for  the 
matrix  material.  A  widely  used  continuum  constitutive  model  for  porous  materials  is  that  of  Gurson 
(1977).  which  has  been  modified  by  Tvergaard  (1982)  with  unit  cell  analysis  to  incorporate  the  ef¬ 
fects  of  void  growth  and  coalescence.  Besides  the  limitations  in  representing  actual  microstructural 
heterogeneities,  a  number  of  these  constitutive  models  do  not  adequately  accommodate  variations 
in  constitutive  parameters  with  evolving  deformation  and  do  not  account  for  post-yield  anisotropy. 
Terada  and  Kikuchi  (1995)  have  tried  to  overcome  this  by  using  the  asymptotic  homogenization  to 
develop  an  extensive  numerical  response  database  in  the  strain  space.  Instantaneous  overall  com¬ 
posite  properties  are  determined  from  discrete  values  of  homogenized  stress-strain  values  at  points 
of  this  database.  This  approach,  however  leads  to  huge  database  to  cover  all  possible  deformation 
paths  and  requires  solving  an  inordinately  large  number  of  RVE  boundary  value  problems.  Fish  et. 
al.  (1997)  have  used  the  idea  of  transformation  strain  fields,  introduced  by  Dvorak  and  Benveniste 
(1992),  to  develop  a  two  point  averaging  scheme  based  on  the  mathematical  homogenization  theory 
with  piecewise  constant  transformation  fields.  However,  approximating  the  eigen-strains  with  low 
order  polynomial  functions  may  not  be  able  to  fully  account  for  large  gradients  in  stresses  and 
strains  between  phases. 

Motivated  by  two  considerations,  a  piecewise  continuous  elastic-plastic  constitutive  model  with 
an  anisotropic  yield  function  is  developed  in  this  work.  The  first  is  an  accuracy  consideration, 
in  that  it  should  account  for  the  microstructural  morphology,  e.g.  spatial  distributions,  shapes, 
sizes  and  properties  of  the  individual  phases,  phase  interactions,  as  well  as  the  evolving  stress  and 
strain  fields.  This  can  be  achieved  if  the  model  is  developed  from  detailed  finite  element  analyses 
of  the  RVE  (e.g.  VCFEM  analysis),  subjected  to  a  wide  variety  of  loading  conditions.  The  second 
is  an  efficiency  consideration,  since  the  creation  of  a  prohibitively  large  numerical  database  with 
a  very  large  number  of  numerical  experiments  is  of  no  consequence.  The  efficient  development 
of  a  constitutive  model,  accounting  for  underlying  evolution  of  state  variables,  is  accomplished  by 
generating  piecewise  continuous  model  peirameters  from  data  in  a  discretized  strain  space  (see  figure 
4.4).  Numerical  data  points  in  the  strain  space  are  systematically  created  through  a  sequence  of 
computational  RVE  analyses  subject  to  an  ordered  set  of  macroscopic  strains  and  strain  paths.  The 
strain  space  in  figure  4.4  is  discretized  into  cubic  elements,  each  containing  32  nodes  or  data  points. 
From  the  computational  RVE  analysis,  constitutive  parameters  like  yield  function  coefficients  and 
plastic  work  are  generated  for  each  nodal  point.  The  constitutive  relation  at  any  point  in  the 
strain  space  are  then  obtained  by  interpolating  nodal  values  using  conventional  shape  functions. 
Elastic-plastic  models  developed  in  the  ensuing  sections  are  for  plane  strain  assumptions. 


Linear  Elasticity 

Orthotropic  homogenized  elastic  material  properties  are  obtained  by  asymptotic  homogenization 
in  conjunction  with  the  VCFE  analysis  of  the  composite  and  porous  microstructui^s  from  equation 
(4.19)  as  explained  in  Ghosh  et.  al.  (1995).  With  plane  strain  assumptions,  three  separate  VCFE 
analyses  are  conducted,  each  corresponding  to  an  independent  component  of  the  macroscopic  strain 
{e„,  Cyy,  exy}.  The  orthotropic  elasticity  tensor  is  stored  for  macroscopic  analysis. 
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4.5.2  Elasto-Plasticity  with  Anisotropic  Yield  Function 

The  inclusion  phase  in  composites  are  assumed  to  be  linear  elastic,  while  the  matrix  phase  is 
assumed  elastic-plastic  for  both  composite  and  porous  materials.  In  plane  strain  modeling,  an 
assumption  that  the  total  plastic  strain  in  the  out  of  plane  or  ‘third’  direction  is  zero,  is  made. 
The  yield  function  can  then  be  described  in  terms  of  the  macroscopic  in-plane  stress  components 
E„„  and  E^v).  The  yield  function  for  porous  and  composite  are  written  using  Hill’s  1948 
anisotropic  yield  function  (see  e.g.  Hill  (1965),  Xie  and  Adams  (1995))  in  conjunction  with  the 
hydrostatic  stress  dependent  models  of  the  Tvergaard-Ourson  (see  e.g.  Gurson  (1977),  Tvergaard 
(1982))  as: 

where  Sxxi^yyi^xy  I'hc  macroscopic  in-plane  strains,  C{6xx1^yy^^xy^^p)  is  a  strain  dependent 
yield  surface  parameter  and  Yf{Wp)  is  the  flow  stress  in  shear.  For  the  composite  materials  the 
dependence  of  pressure  on  yielding  is  deemed  negligible  and  the  hydrostatic  stress  coefficient  H 
is  ignored.  Coefficients  C{exx.eyy,exy,Wp),  Yf{Wp)  and  H  are  determined  from  computadonal 
experiments  detailed  next.  The  increment  of  plastic  strain  is  obtained  from  the  yield  function  $ 
by  using  the  associated  flow  rule  for  hardening  materials,  i.e. 

Coefficient  Evaluation 
A.  H  and  YfjWp) 

Computational  exercises  indicate  that  the  variation  of  H  with  increasing  hydrostatic  loading  is 
not  significant.  It  is  therefore  assumed  to  be  a  constant  for  all  load  histories.  This  assumption  is 
consistent  with  the  Tverggard-Gurson  models,  where  H  is  determined  in  terms  of  the  initial  void 
volume  fraction.  The  constitutive  parameters  H  and  Yf{Wp)  are  evaluated  in  a  coupled  manner  by 
solving  the  microstructural  RVE  boundary  value  problems  with  two  distinct  loading  conditions  viz. 
(i)  biaxial  tension  loading  (Sn  =  Eyy  =  Ejy  =  0)  and  (ii)  pure  shear  loading  (Eiy  =  E^/,, 

Eli  =  Syy  =  0)-  For  load  condition  (i),  equation  (4.26)  becomes: 

^i^kyd,  ^hyd,  0,  Wp)  =  Hcosh  -1  =  0  (4.27) 


and  for  load  condition  (ii),  it  becomes: 

4.(0,0,E,/„H'p)  =  jp|y  +  H-l=0  or  (‘‘•28) 

The  values  of  H  and  Yf{Wp)  are  determined  iteratively  from  equation  (4.28)  and  further  validated 
against  equation  (4.27).  The  steps  are  as  follows. 

1.  Solve  a  macro- micro  boundary  value  problem  with  RVE  homogenization,  with  incremental 
pure  shear  loading.  Obtain  macroscopic  plastic  work  by  averaging  the  the  microstructural 
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plastic  work  (  Wp  =  /q  crijifjdQ  )  and  plot  the  macroscopic  shear  stress  as  a  function  of 
plastic  work  Wp. 

2.  Assume  a  starting  value  for  H  (e.g.  3*/o  as  in  Gurson  (1977),  Tvergaard  (1982)  and  evaluate 
Yf{Wp)  from  equation  (4.28). 

3.  Solve  a  pure  macroscopic  boundary  value  problem  with  incremental  biaxial  loading,  using  the 
homogenized  elastic-plastic  constitutive  relation  and  associated  flow  rule  with  yield  function 
(4.26). 

4.  Plot  the  Cix  -  S/iyd  and  Cyy  -  E/jyd  curves  for  the  entire  history  of  biaxial  loading. 

5.  Solve  a  macro-micro  boundary  value  problem  with  RVE  homogenization  with  the  same  in¬ 
cremental  biaxial  loading  as  in  the  previous  step.  Plot  the  Hhyd  -  ^xx  and  T^hyd  -  ^yy  curves 
for  the  entire  history  of  loading. 

6.  Compare  the  results  of  steps  4  and  5.  If  the  two  curves  from  both  methods  are  within  a  preset 
tolerance  everywhere,  then  the  process  is  terminated  and  value  of  H  in  step  2  is  accepted. 
Otherwise,  the  entire  sequence  is  repeated  with  a  different  value  of  H. 


B.  C{e,j,Wp) 


The  coefficient  C  in  equation  (4.26)  is  found  to  vary  considerably  with  evolution  of  plastic  de¬ 
formation  and  examples  of  its  variation  with  straining  and  plastic  work  are  shown  in  figures  4.5. 
While  it  is  assumed  to  be  a  function  of  the  total  strain  and  plastic  work,  its  dependence  on  load 
history  is  assumed  to  be  negligible.  In  the  discretized  strain  space  of  figure  4.4,  the  value  of  a  piece- 
wise  continuous  C  at  any  point  may  then  be  obtained  by  interpolation  from  nodal  values  according 
to 

32 

C7(eii,  Cyy,  Ciy,  Wp)  =  Ca{Wp) Nai^xx  1  ^yy  1  ^xy)  (4.29) 

a=l 

where  Cq  are  the  nodal  values  and  Na  are  shape  functions  for  a  32-noded  brick  element. 
Generation  of  Discretized  Strain  Space  and  Nodal  Parameters 

The  nodal  values  of  macroscopic  stresses  (Sn,  Eyy,  Exy)  and  the  corresponding  plastic  work  Wp  are 
first  evaluated  at  each  nodal  point  in  a  subspace  of  the  Cxi  —  Syy  —  Ciy  space  by  solving  incremental 
macro-microscopic  boundary  value  problems  with  VCFEM  and  asymptotic  homogenization.  In 
this  process,  macroscopic  strain  increments  are  applied  to  the  RVE  subjected  to  periodic  bound¬ 
ary  conditions  (see  e.g.  Ghosh  et.  al.  (1996,1997)).  Strain  increments  are  applied  along  a  radial 
line  in  the  strain  space,  such  that  a  constant  ratio  between  strain  components  is  maintained,  i.e. 
Acii  :  Acyy  :  Aciy  =  1  :  tan6  :  (1  +  tav?9)tan(l>,  where  9  and  (p  are  the  angulac  coordinates  in 
the  strain  space  of  figure  4.4.  The  fiow  stress  Yf(Wp)  at  each  node  in  figure  4.4  can  be  obtained 
from  the  shear  stress-plastic  work  plot  and  equation  (4.28).  Prom  the  values  of  liiacroscopic  vari¬ 
ables  (Sxi,Syy,Sjy,y/(Wp))  at  a  node  corresponding  to  the  end  of  an  increment,  the  coefficient 
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C(eii,ej,j„eiy,  Wp)  is  calculated  using  equation  (4.26). 

From  the  symmetry  conditions,  only  a  quarter  of  the  Cn  —  e.yy  —  exy  strain  space  is  considered 
for  loading  such  that  0°  <  9  <  180"  in  the  Cu  -  Cyy  plane  and  0"  <  0  <  90"  outside  of  this  plane. 
This  is  indicated  in  figure  4.4.  This  chosen  subspace  of  the  strain  space  is  divided  into  16  cubic 
brick  element  with  32  nodal  points  each.  The  location  of  elements  in  the  strain  space  are  selected 
to  optimally  account  for  the  variations  in  C.  These  variations  in  C  with  the  coordinate  angle 
9  (location  in  the  Cxx  -  eyy  plane)  and  the  plastic  work  Wp  for  the  different  microstructures  are 
plotted  in  figure  4.5.  The  parabolic  form  of  C  in  figure  4.5a  is  consistent  with  the  quadratic  term 
(Eix  -  Syy)^  in  the  yield  function.  The  minimum  values  occur  near  9  =  135"  corresponding  to  a 
pure  deviatoric  state.  The  coefficient  subsequently  increases  to  account  for  the  increase  in  plastic 
work  in  the  yield  function  In  figure  4.5b,  the  coefficient  C  as  a  function  of  the  plastic  work, 
which  corresponds  to  the  radial  direction  in  the  strain  space,  is  plotted.  The  value  of  C  stabilizes 
beyond  a  value  of  the  plastic  work,  which  is  used  as  the  outer  boundary  of  the  strain  space  envelope 
in  figure  4.4. 


4.5.3  Numerical  Implementation  of  the  Constitutive  Model 

The  elastic-plastic  constitutive  model  for  composite  and  porous  materials  is  derived  from  the 
anisotropic  yield  function  (4.26)  with  associated  flow  rule  and  isotropic  hardening.  In  an  incre¬ 
mental  form,  the  stress  increments  AE^  are  related  to  elastic  increments  of  strains  (Ae*;  -  Ae^,) 
admitting  additive  decomposition,  as 

AEy  =  F;^,/(Ae,,  -  Ae^,)  (4.30) 

where  is  the  homogenized  elasticity  tensor.  Using  associated  flow  rule,  components  of  the 
plastic  strain  increment  are  obtained  as: 


(4.31) 


Elimination  of  the  flow  parameter  AA  from  the  above  equations  results  in  the  two  equations 


-  Ae^ 


-  Ae^ 


(4.32) 


These  equations  are  solved  using  the  backward  Euler  integration  method,  with  gradients  evaluated 
at  the  end  of  the  increment.  With  known  increments  of  strain,  the  resulting  set  of  equations  (4.32) 
together  with  the  yield  function  (4.26)  are  solved  iteratively  by  using  the  Newton-Raphson  method. 
The  stress  increments  are  obtjiined  by  the  following  steps. 

1.  Initialize  values  of  AEn,  AEyy  arid  AEiy.  c. 


2.  Calculate  the  gradient  (^)  of  the  yield  function  and  solve  for  the  increments  of  plastic 
strain  Ae^j,  ^®yy  smd  Ae^y  from  equations  (4.30)  and  (4.32).  Update  the  stresses  and 
plastic  work  using  the  relation  AWp  =  EnAe^j.  -t-  Eyy  ^^yy  ^ly^^iy 
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3.  If  $  <  toll  and  correction  to  plastic  strain  increment  5^^  <  tol2,  where  toll  and  tol2  are 
prescribed  tolerances,  then  stop.  Otherwise  go  to  step  2. 

The  parameter-C  is  then  obtained  from  the  interpolation  equation  (4.29). 

4,5.4  Numerical  Examples  with  the  Continuum  Constitutive  Model 

The  elastic-plastic  constitutive  model  Is  validated  by  comparing  the  results  of  macroscopic  numerical 
simulaiions  with  those  generated  by  macro-micro  scale,  analysis  using  asymptotic  homogenization. 
Examples  are  conducted  for  both  composite  and  porous  materials  with  different  microstructure 
niQrphologies,  viz.  different  shapes,  sizes  and  spatial  distributions  of  the  heterogeneities. 

Analysis  of  Composite  Microstructures 

Six  microstructural  RVE’s  of  20%  V/  Alumina-aluminum  composite  with  Alumina  fiber  in  alu¬ 
minum  matrix,  as  shown  in  figure  4.6,  are  analyzed.  The  RVE’s  are  classified  as: 
laj:  Square  edge  pattern  with  a  circular  inclusion  (Cl) 

(bj:  Square  edge  pattern  with  an  elliptical  inclusion  (aspect  ratio  |  =  3)  (C2) 

(c) :  Random  pattern  with  25  identical  circular  inclusions  (C3) 

(d) :  Horizontally  aligned  random  pattern  with  25  identical  elliptical  inclusions  (C4) 

(e) :  Randomly  oriented  random  pattern  with  25  identical  elliptical  inclusions  (C5) 

(f) :  Random  pattern  with  17  random  shape  and  size  inclusions  (C6) 

The  material  properties  for  the  elastic  Alumina  fiber  are: 

Young’s  Modulus  (Ec)=  344.5GPa,  Poisson  Ratio  (i/c)=  0.26; 
and  for  the  elastic-plastic  Aluminum  matrix  are: 

Young's  Modulus  {^Erri)—  68.3  GPa,  Poisson  Ratio  (i/',7j)=0.30.  Initial  Yield  Stress  (^0).  55  MPa, 
and  Post  Yield  hardening  law:  ^eqv  —  Yo  +  2.08eP'„. 

The  RVE's  are  subjected  to  four  different  types  of  loading  viz. 

•  LI:  Pure  shear  loading  with  increments ( A En  =  AEyy  =  0.  AE^y  =  AE) 

•  L2:  Uniaxial  tension  loading  with  increments  (AEji  =  AE.  AEyy  =  AEiy  =  0) 

•  L3:  Bisucial  tension  loading  with  increments  (AExx  =  =  AE,  AExy  =  0) 

•  L4:  Biaxial  tension-compression  loading  with  increments  (AEn  =  —AEyy  =  AE,  AExy  =  0) 

The  stress-strain  response  of  the  six  composite  microstructural  RVE’s  with  the  four  loading 
conditions  are  conducted  and  the  results  for  simple  tension  (LI)  are  depicted  in  figure  4.7.  The 
results  by  the  constitutive  model  and  two-scale  asymptotic  homogenization  approach  are  generally 
found  to  agree  very  well  for  the  entire  range  of  loading  upto  fairly  high  level  of  straining.  The 
only  discrepancy  is  found  with  the  biaxial  tension  loading  condition  (L3),  for  which  the  devia¬ 
tion  strains  sire  shown  in  table  1.  However  the  deviations  occurs  at  high  strains  levels,  for  which 
the  stresses  are  nearly  twice  the  matrix  yield  stress.  It  is  important  to  note  that  the  deviation 
of  continuum  model  response  from  the  two-scale  asymptotic  homogenization  response  can  be  used 
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as  a  signal  for  switching  from  the  fornisr  to  the  latter  type  of  analysis  in  a  multiple  scale  simulation. 

As  an  example  of  a  structural  analysis  with  the  two  approaches,  a  square  plate  with  a  square 
hole  is  solved  with  tension  loading  on  two  opposite  faces.  A  quarter  of  the  plate  with  symmetry 
and  loading  conditions  is  shown  in  figure  4.8a.  A  total  traction  of  55  MPa  is  applied  in  10  equal 
increments.  The  material  is  a  20  %  V/  Alumina-aluminum  composite  with  the  microstructural 
RVE  consisting  of  15  identical  circular  inclusions  dispersed  randomly  in  the  matrix  (figure  4.8b). 
The  material  properties  are  the  same  as  in  the  previous  example.  Comparison  of  the  results  are 
made  through  contour  plots  of  the  macroscopic  stress  Exx  (shown  in  figure  4.9)  and  macroscopic 
plastic  work  (not  shown).  The  figures  reveal  that  at  all  locations  the  difference  in  the  two 
approaches  is  less  than  1%.  However,  the  computational  efficiency  of  the  macroscopic  analysis 
with  the  continuum  constitutive  model  is  far  superior  than  the  two  scale  analysis.  The  scale  up 
in  efficiency  for  this  problem  is  approximately  75000%,  and  is  therefore  very  desirable  when  only 
macroscopic  results  are  of  interest. 

Analysis  of  Porous  Microstructures 

The  six  microstructural  RVEs  are  analyzed  again  for  porous  materials,  with  voids  replacing  the  in¬ 
clusions  in  figure  4.6.  The  material  considered  is  an  aluminum  alloy  with  20%  void  volume  fraction 
and  the  following  elastic-plastic  properties; 

Young’s  Modulus  (£'m)=  68-3  GPa,  Poisson  Ratio  {um)=  0.30,  Initial  Yield  Stress  (yo)=  55  MPa, 
and  Post  Yield  hardening  law:  ae,v  =  Vo  +  The  same  four  load  histories  (L1,L2,L3,L4) 

are  applied.  An  important  difference  between  the  composite  and  porous  microstructures,  is  that  in 
the  latter  plastic  strain  localization  in  small  regions  is  a  common  occurrence  depending  on  the  void 
morphology  and  the  nature  of  loading.  Such  nonhomogeneous  distribution  of  plastic  strains  is  a 
major  source  of  discrepancies  between  responses  by  the  two  approaches  and  act  as  limiters  for  the 
range  of  application  of  the  continuum  model.  Microstructures  with  homogeneous  and  nonhomoge¬ 
neous  distributions  of  plastic  strain  are  shown  in  figure  4.10.  For  the  microstructure  VI  (square 
edge  pattern  with  a  circular  void)  the  strain  distribution  is  quite  uniform  in  pure  shear  loading, 
while  for  the  microstructure  V2  (square  edge  with  an  elliptical  void)  there  is  intense  localization 
with  narrow  ligaments.  Consequently  the  continuum  model  ceases  to  be  effective. 

As  in  the  case  with  composites,  the  main  challenge  for  the  homogenized  constitutive  model 
is  encountered  during  simulations  with  biaxial  tension  loading,  i.e.  (AEn  =  AEyy  =  AEj,  and 
AEiy  =  0).  The  first  term  in  the  yield  function  (4.26)  drops  out  for  this  loading  and  the  model 
delivers  the  same  amounts  of  plastic  strains  in  the  i  and  y  directions.  Due  to  the  lack  of  anisotropy 
in  the  hydrostatic  term  in  the  yield  function,  the  continuum  model  is  effective  only  for  those 
microstructimes  that  exhibit  near-isotropic  plastic  behavior  for  this  loading.  Microstructures  V2 
(square  edge  with  elliptical  void)  and  V4  (horizontally  aligned  elliptical  voids  with  random  spatial 
distribution)  shows  very  different  strain  responses  for  each  direction  with  biaxial  loading.  Thus  the 
continuum  constitutive  model  is  largely  ineffective  for  these  RVE’s.  The  microstructures  V5  (ran¬ 
domly  oriented  identical  elliptical  voids  with  random  spatial  distribution)  and  V6  (random  spatial 
distribution  with  random  shape  and  size)  also  show  significant  plastically  induced  directional  effects 
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and  the  constitutive  models  are  therefore  restricted  to  the  elastic  range.  The  list  of  performance 
and  strain  ranges  of  all  the  microstructures  with  the  different  loading  are  given  in  table  2.  The 
microstructures  VI  and  V3  exhibit  relatively  isotropic  responses  and  yield  satisfactory  agreement 
between  responses  by  the  constitutive  model  and  the  two-scale  asymptotic  homogenization.  Com¬ 
parisons  of  stress-strain  responses  for  biaxially  and  uniaxially  loaded  RVE’s  are  made  in  figures 
4.11.  These  show  very  good  agreement.  The  RVE’s  V2  and  V4  exhibit  intense  localization  early 
on  in  the  straining,  while  the  RVE’s  V5  and  V6  exhibit  marked  anisotropy  with  biaxial  loading. 
Plastic  flow  predictions  for  these  RVE’s  with  the  continimm  model  are  therefore  not  reliable. 

4.5.5  Level-0/ 1  Mesh  Enrichment  by  h- Adaptation 

The  transition  between  various  levels  in  this  model  is  augmented  by  an  adaptive  mesh  refinement 
involving  element  subdivision  or  h-refinement  in  the  level-0  and  level-1  regions.  This  local  mesh 
enrichment  is  intended  to  serve  two  purposes,  viz.  (a)  to  identify  and  reduce  discretization  error 
in  the  computational  model  and  (b)  to  reduce  modeling  error  by  zooming  in  on  regions  of  evolving 
localization  due  to  microscopic  non-homogeneity.  The  latter  is  also  effective  in  bridging  the  gap 
between  the  macro-  and  micro-  scale  elements  by  successive  element  refinement  in  critical  regions, 
as  shown  in  figures  (4.2a  and  4.16).  The  adaptive  h-refinement  has  been  extensively  discussed  in 
literature  (see  e.g.  Zhu  and  Zienkiewicz  (1988,1992),  Melosh  and  Marcal  (1977),  Demkowicz  et.  al. 
(1985),  Bass  and  Oden  (1987))  and  the  latest  advances  in  adaptive  methods  in  mechanics  have  been 
presented  in  a  compilation  by  Ladev‘eze  and  Oden  (1998).  Within  the  category  of  h-adaptive  pro¬ 
cedures,  two  methods,  viz.  the  mesh  enrichment  or  refinement  methods  (e.g.  Zhu  and  Zienkiewicz 
(1988),  Bass  and  Oden  (1987))  and  mesh  regeneration  methods  (e.g.  Zhu,  Hinton  and  Zienkiewicz 
(1993),  Paulino  et.  al.  (1999))  have  evolved.  While  the  mesh  regeneration  methods  have  been 
preferred  for  their  relatively  higher  efficiency,  the  enrichment  method  is  deemed  more  suitable  for 
the  present  work.  This  is  due  to  the  fact  that  the  regeneration  method  alters  element  locations 
with  respect  to  the  underlying  material  RVE  and  will  necessitate  cumbersome  mapping. 

Following  procedures  outlined  in  Bass  and  Oden  (1987),  the  mesh  refinement  procedure  entails 
subdividing  each  quadrilateral  QUAD4  element  into  four  smaller  elements,  thus  adding  four  new 
nodes  on  the  boundairy  zmd  one  in  the  center.  For  compatibility  with  the  displacement  fields  on  the 
boundary  of  larger  adjacent  elements,  linear  constraint  relations  are  imposed  on  the  new  boundary 
nodes  of  the  subdivided  element.  A  static  condensation  process  is  then  used  to  eliminate  the  explicit 
inclusion  of  the  new  nodes  in  the  global  stiffness  and  load  vectors  (see  Bass  and  Oden  (1987)  for 
details).  Numerical  integration  in  each  QUAD4  element  after  adaptation  is  performed  by  one  point 
reduced  integration  with  hourglass  control  (Koh  and  Kikuchi  (1987)). 

4.6  Computational  Subdomain  Level-1  in  the  Hierarchical  Model 

The  level-1  regions  (fla  C  ^mat{p))  the  computational  domain  cire  intended  as  ‘swing’  or  ‘tran¬ 
sition’  regions,  where  microscopic  information  in  the  RVE  is  used  to  decide  whether  microscopic 
computations  are  necessary  for  these  regions.  They  are  identified  by  locally  high  gradients  of  macro- 


65 


scopic  variables  e.g.  stresses,  strains  strain  energy  etc.  that  are  indicators  of  imminent  damage 
evolution  or  localization  in  the  microstructure.  Computations  in  these  regions  cire  still  based  on 
assumptions  of  macroscopic  uniformity  and  periodicity  of  the  RVE.  Concurrent  with  macroscopic 
simulations,  computations  are  necessary  in  the  the  microstructure  to  monitor  the  initiation  and 
growth  of  damage.  The  RVE  response  for  level-1  elements  is  explicitly  calculated  using  asymptotic 
homogenization  and  VCFEM  analysis.  The  computations  in  elements  belonging  to  the  level-1  do¬ 
main  undergo  a  sequence  of  finite  element  analyses  as  follows. 

(a)  Microstructural  analysis  of  the  RVE,  subjected  to  the  sequence  of  unit  macroscopic  strains  or 
increments  with  periodicity  boundary  conditions  to  generate  homogenized  tangent  modulus 

(b)  Macroscopic  analysis  of  the  structure  using  of  step  (a),  to  evaluate  macroscopic  vari¬ 

ables  e.g.  stresses  and  strains  due  to  applied  loads. 

(c)  Microstructural  analysis  of  the  RVE  at  each  sampling  point  (e.g.  integration  point)  of  macroscopic 
elements,  with  actual  macroscopic  strains  and  increments  obtained  from  step  (b)  and  periodicity 
conditions  on  the  RVE  boundary.  Microscopic  stresses  and  strains  (<T^€‘)  Jire  thus  calculated  in 
the  RVE's  of  each  element.  For  linear  elastic  constituent  phases,  the  microstructural  stress  recovery 
process  can  be  achieved  by  using  a  linear  combination  of  solutions  with  unit  strains  from  step  (a). 
However,  explicit  solution  of  nonlinear  equations  are  required  for  problems  with  nonlinearity.  This 
is  executed  in  an  iterative  manner  in  this  work. 

Computational  requirements  of  elements  in  this  level  are  considerably  higher  than  that  for  level 
0.  At  each  integration  point  in  the  macroscopic  computational  mesh,  a  complete  microstructuraJ 
analysis  of  the  RVE  problem  is  done  several  times  (at  least  3  times  for  step  (a)  and  once  for  step  (c) 
in  the  present  case)  within  each  iteration  loop.  Thus  the  level- 1  elements  are  computationally  more 
expensive  compared  to  level-0  elements.  It  is  therefore  important  to  design  robust  criteria  to  avoid 
redundant  element  transition  from  level-0  to  level- 1.  It  is  also  critical  that  this  level  be  discontinued 
once  the  microstructural  damage  or  instability  evolved  beyond  pre-determined  threshold  values. 

4.7  Computational  Subdomain  Level-2  in  the  Hierarchical  Model 

Level-2  regions  (figure  4.2b)  axe  classified  as  those  with  severe  microstructural  nonuniformities  in 
the  form  of  evolving  damage.  This  results  in  loss  of  statistical  periodicity  of  the  assumed  RVE  and 
these  regions  may  be  identified  with  ^MAT[np)  equation  (4.5).  In  the  computational  model,  the 
level-2  elements  (P/j  C  ^MAT{np))  materialize  from  the  microstructure  of  level-1  elements.  It  is  as¬ 
sumed  that  prior  to  this  transition,  the  level- 1  elements  have  been  refined  to  reach  sufficiently  high 
spatial  resolution.  In  0(21  the  macro-micro  model  of  level- 1  switches  to  a  completely  microscopic 
model. 

The  method  of  generating  the  level-2  microstructure  in  each  element  is  illustrated  in  figure 
(4.12).  An  extended  microstructure  is  first  created  by  repeating  RVE’s  in  succession,  to  cover  the 
entire  region  of  the  macroscopic  level- 1  elements  in  transition  to  level-2.  In  equation  (4.7),  a  local 
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non-periodic  region  is  created  as  : 


yMATinp)  ^  (4  33) 

where  corresponds  to  the  RVE’s  in  a  periodic  domain,  adjusted  for  microscopic  evolution.  This 
is  then  overlayed  by  the  macroscopic  elements  to  accurately  encase  the  level-2  region  with  clearly 
delineated  boundaries.  Each  level-2  element  now  contains  a  heterogeneous  material  distribution 
(Y^)  that  is  defined  as  the  intersection  of  the  non-periodic  material  region  and  the 

element  domain  Q^2 

ye  ^  y,MAr(np)  ^  (4  34) 

This  region  is  subsequently  tessellated  into  a  mesh  of  Voronoi  cell  elements  (figure  4.12b).  Trac¬ 
tion  continuity  between  level-2  microstructural  region  and  neighboring  level-0/ level- 1  elements  is 
incorporated  through  a  layer  of  transition  elements. 


4.7.1  Damage  Criterion  for  Particle  Cracking  in  Level-2 

The  level-2  VCFEM  modeling  consist  of  brittle  reinforcing  particles  and  a  ductile  matrix  material. 
For  the  brittle  particulate  materials,  microstructural  damage  initiation  is  assumed  to  be  governed 
by  a  maximum  principal  stress  based  criterion.  In  this  criterion,  a  crack  is  initiated  when  the 
maximum  principal  stress  in  tension  exceeds  a  critical  fracture  stress  at  a  point.  In  the  compu¬ 
tational  procedure,  complete  particle  splitting  is  assumed  to  occur  in  the  form  of  an  elliptical  void, 
normal  to  the  principal  direction,  as  soon  as  the  principal  tensile  stress  reaches  In  the  case 
of  particle  splitting,  the  crack  tip  extends  nominally  into  the  matrix.  In  the  incremental  compu¬ 
tational  procedure,  more  than  one  point  may  exceed  the  critical  cr"’  value  during  increment.  The 
location  of  a  single  crack  is  determined  by  a  weighted  averaging  method  as: 


^damage 


^  c^x,y) 


Vdamagt  — 


Ey 


Cycr 


(4.35) 


where  cr;(i,  y)  corresponds  to  all  values  of  maximum  tensile  principal  stress  larger  than  o-”"  in  the 
particle.  The  crack  is  oriented  at  right  angles  to  the  principal  stress  directions  at  (xdamagei  Vdamage) 
and  extends  to  the  interface  on  both  sides. 


4.8  Coupling  the  Levels  in  the  Hierarchical  FE  Model 

While  level-0  elements  {Eiq  €  fiio)  stud  level- 1  elements  {Ei\  €  flu)  are  coupled  naturally  through 
the  familiar  assembly  process,  the  interf2icing  of  level-2  [Ea  €  fl(2)  elements  with  either  of  the  first 
two  requires  more  attention.  The  mismatch  in  the  number  of  boundary  nodes  in  these  elements 
necessitate  the  introduction  of  transition  elements  (E(r  €  fl/a).  acting  as  buffer  zones  as  shown  in 
figure  4.13.  Both  Ei2  and  Etr  elements  employ  VCFEM  for  setting  up  the  element  Stiffness  matrices 
and  load  vectors.  The  entire  computational  domain  is  thus  comprised  of 

n'  =  {n,o  U  fin  u  Qi2  :  (iiQ  =  ^kL\Eio-,  fill  =  Q,2  =  ^kL\Ei2  u 
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for  which  the  nodes  are  differentiated  as  (see  figure  4.13): 

(i)  (nd*^^)  or  internal  level-2  Voronoi  cell  element  nodes  that  are  not  on  any  interface  or  boundary, 

(ii)  (nd^^''^)  or  Voronoi  cell  element  nodes  on  the  E12-E12  interface, 

(iii)  or  Voronoi  cell  element  nodes  on  Ei2-Etr  interface, 

(iv)  (nd'®^/*’’)  or  nodes  on  Eio/i-Etr  interface  that  belong  to  Eio  and  En  elements, 

(v)  (nd*‘^)  or  Voronoi  cell  element  nodes  on  E/o/r^^tr  interface  that  do  not  belong  to  Eio  and  En 
and  need  to  be  statically  condensed.  In  an  incremental  analysis  for  elasto-plasticity,  the  principle 
of  virtual  work  for  the  computational  domain  at  the  end  of  the  n  -  th  increment  may  be  written 
as  a  scalar  valued  function  G  in  terms  of  variables  at  different  levels  as: 


G”+^(Au,(iu)  = 
-H 

-f 


tiSuidT  -  f  tiSujdr  -  [  tiSuidT 

Jr?*'  Jr?*' 


(4.36) 


In  an  iterative  solution  process  with  the  Newton-Raphson  method,  a  consistent  linearization  by 
taking  the  directional  derivative  of  along  incremental  displacement  vector  Au,  yields  the 

tangent  stiffness  matrix.  For  the  i  -  th  iteration,  the  assembled  equations  have  the  following 


structure. 

r  ;^/o/i./o/i  j^io/u2yj  ^ijion  )'_(  1 

^(2,/0/l  I  ^ul2B  I  -  I  ^pl2  I 


Here  corresponds  to  displacements  at  nodes  (nd'°/')  that  belong  to  elements  Eio  and  En 

in  the  computational  region  U  fiji  as  shown  in  figure  4.13.  It  should  be  noted  that  they  also 
include  nodes  at  the  interface  of  elements  EiojEn  and  elements  Etv  The  displacements  AC/'^®  on 
the  other  hand  corresponds  to  nodes  (nd'^^‘'')  on  the  interfaces  of  elements  Etr  and  elements  En  or 
to  nodes  (nd‘^^^)  on  the  interfaces  between  two  En  elements.  Contributions  to  the  stiffness  matrix 
[F]  and  load  vector  {F}  for  elements  in  U  flij  may  be  obtained  as 


/0/l,/0/ 


1)«=  f 
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f(Un  +  Au‘)- 


a,0  correspond  to  the  node  numbers  and  Nq  are  the  shape  functions  of  macroscopic  elements. 
For  the  elements  Etr  and  En  belonging  to  Cl^,  contributions  to  the  stiffness  matrix  and  load 
vectors  are  obtained  by  VCFEM  calculations,  together  with  static  condensation.  The  transition 
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element  facilitates  continuous  variation  from  microscopic  to  macroscopic  elements  without  jumps 
or  discontinuities.  On  the  £'/o/i-£'tr  interfaces,  this  is  accomplished  by  constraining  displacements 
at  nodes  to  conform  to  displacement  interpolation  of  the  adjacent  Eiq  or  En  elements.  The 
constraints  at-the  nodes  are  applied  in  terms  of  displacements  at  as 

{AC/*=}  =  (4.38) 

where  [Q]  is  the  constraining  matrix.  For  bilinear  QUAD4  level-0/1  elements,  each  row  of  [Q] 
consists  of  the  inverse  of  the  distance  of  the  constrained  node  to  the  corner  nodes.  The  interfaces 
with  the  Ei2  elements,  i.e.  ihb  En-Etr  interfaces  are  treated  in  the  same  way  as  E12-E12  interfaces. 


The  displacements  AC/'^^  in  equation  (4.37)  correspond  to  those  at  the  boundary  nodes  (nd'*^^. 
nd'"'"’,  nd^°'''‘'’  and  nd**^)  of  level-2  elements  that  contain  the  microstructural  VCFE  model.  To 
account  for  the  contribution  of  the  internal  nodes  (nd'^^),  it  is  therefore  necessary  to  use  static 
condensation  and  recovery  process  for  representing  the  VCFEM  displacement  solutions  at  the 
internal  nodes  in  terms  of  the  displacements  at  boundary  nodes  AC/'^®,  where 
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(4.39) 


where  [I]  and  [Q]  are  the  identity  and  constraint  matrices  respectively.  In  a  VCFEM  solution 
process  the  stiffness  matrix  and  the  load  vector  may  be  partitioned  accordingly  as 


•  f^l2B,l2I  I  j  ^^128  1  J  ^pl2B  1 

K^2I,12B  J^12I,12I  J  I  ^^l2I  I  -  I  ^fl2I  | 


(4.40) 


Static  condensation  of  the  internal  degrees  of  freedom  yields 


(4.41) 

These  stiffness  matrices  and  load  vectors  axe  then  used  in  the  assembly  process  of  equation  (4.37). 


4.9  Adaptation  Criteria  for  Various  Levels 

It  is  evident  that  appropriate  criteria  need  to  be  established  for  transitioning  the  computational 
subdomains  from  one  level  to  another.  In  addition  to  level  transitions,  element  refinement  by  h- 
adaptation  is  also  executed  in  level-0  and  level- 1  regions  for  reducing  discretization  CTror,  identifying 
regions  of  high  solution  graidients  and  zooming  in  to  reduce  the  scale  gap  between  macroscopic  and 
microscopic  elements.  The  criteria  used  for  element  refinement  and  level  transition  axe  delineated 
below. 
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•  Level  0/1  h-adaptation  may  be  executed  based  on  an  error  measure  may  be  defined  as 

^,  =  i|/(u-Uh)||  (4.42) 

where  u  is  the  true  solution,  Uh  is  its  finite  element  approximation,  /  is  any  appropriate  function 
of  deformation  e.g.  energy,  stresses,  strains  etc.,  and  ||  •  ||  denotes  a  norm  of  the  function.  Work 
on  h-method  of  mesh  refinement  (see  e.g.  Zhu  and  Zienkiewicz  (1988).  Melosh  and  Marcal  (1977), 
Demkowicz  et.  al.  (1985))  have  proposed  various  forms  of  /  and  norms  to  optimize  the  compu¬ 
tational  effort  and  improve  reliability.  In  this  work,  elements  are  subdivided  based  on  an  element 
le%-el  parameter  defined  as 

.Il/(uh)ll  .  ,4  ox 

{Ee)cr  =  qi  r^-  (-i^S) 


sfWE 

where  {Ee)cT  corresponds  to  a  threshold  element  parameter,  above  which  the  element  needs  to  be 
refined  and  qi  is  a  prescribed  quality  index.  ||/(uh)||  =  ll/(uh)lle  is  the  sum  of  element  norms 
of  any  function  for  the  entire  computational  domain  and  NE  is  the  number  of  elements.  Such  error 
criteria  is  conventionally  used  to  equi-distribute  errors  due  to  discretization  by  mesh  refinement  for 
relatively  well  behaved  problems.  However  for  mesh  dependent  problems  with  inherent  localizations 
or  singularities  in  the  solutions,  the  above  criteria  may  be  used  for  signaling  steep  gradients  in  the 
solution  variables  and  may  be  used  for  increasing  the  resolution  in  the  region  of  high  local  gradients. 
In  this  work  ||/(uh)|le  is  taken  to  be  the  maximum  difference  in  plastic  work  with  neighboring 
elements,  i.e.: 

l|/(Uh)||e  =  Max\(  [  dfi)^  -  (  [ 

j  Q  •/  n 

with  a  value  of  91=1.5  . 


•  Level-0  to  level- 1  transition  of  Eio  elements  takes  place  in  accordance  with  criteria  based 
on  macroscopic  variables  e.g.  (S,  e)  in  the  continuum  model  that  are  related  to  critical  micro¬ 
scopic  variables.  Additionally,  this  transition  is  activated  when  the  continuum  level  model  in  section 
4.5.1  starts  to  digresses  considerably  from  the  results  of  two-scale  homogenization.  Different  criteria 
are  used  for  composite  and  porous  materials. 

A.  For  composite  microstructure  with  inclusions: 

(E/)  >  E^  =  r*cr"^  and/or  +  ^yy)  ^  (444) 

Here  E/  is  the  maximum  principal  stress  in  Eto,  E'^’’  is  a  critical  stress  that  is  a  pre-determined 
fraction  r  of  the  critical  fracture  stress  [a^^)  for  microstructural  inclusions.  The  fraction  r  is 
selected  to  be  5  in  this  study.  The  second  condition  is  established  since  the  earliest  digression  from 
the  homogenization  results  is  observed  for  biaxial  loading,  and  E^  established  from  the  results  of 

section  4.5.4.  ^ 

B.  For  porous  microstructure  with  voids:  Strain  based  criteria  are  deemed  more  important  in  the 
case  of  damage  in  porous  materials  and  hence  transition  is  activated  when 

(eP)>(eP)"=r*(6P)^  and/or  (4.45) 
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where  €’’{=  is  the  macroscopic  effective  plastic  strain  in  Eio  and  {e^}^  is  a  pre-determined 

fraction  r  of  a  suitably  assessed  high  microstructural  plastic  strain  (eP)".  Again  is  established 
from  the  limiting  values  of  biaxial  strain  in  section  4.5.4. 

•  Level-1  to  leveI-2  transition  of  En  elements  occur  when  a  sufficiently  high  spatial  resolu¬ 
tion  is  attained  by  h-refinement  and  when  the  RVE  is  deemed  to  be  on  the  verge  of  significant 
damage  evolution.  The  adaptive  criteria,  which  monitor  the  transition  from  elements  En  to  E12 
are  prescribed  in  terms  of  topological  evolution  of  microscopic  damage  as: 

A.  For  composite  microstructure  with  inclusions: 


Rdamage  — 

>  R 

(4.46) 

B.  For  porous  microstructure  with  voids; 

ADR 

(4.47) 

Rdamage  ~ 

AR 

where  Rdamage  is  the  ratio  of  the  number  of  damaged  inclusions  (NPD)  to  the  total  number  of 
inclusions  (NP)  for  composites.  For  porous  microstructures,  it  is  the  ratio  of  the  damaged  area 
corresponding  to  regions  which  have  high  plastic  strains  to  the  total  RVE  area.  i?cr  is  a  prescribed 
critical  ratio  and  varies  from  problem  to  problem. 

4.10  Numerical  Examples  with  the  Adaptive  Multilevel  Model 

Numerical  examples  are  solved  to  understand  the  effect  of  structural  geometry  and  microstructural 
morphologies  on  the  macroscopic  and  microscopic  responses.  In  all  examples,  the  inclusions  are 
assumed  to  be  linear  elastic  which  can  crack  by  the  principal  stress  criterion  and  the  matrix  is  elastic- 
plastic.  In  the  first  example,  a  RVE  consisting  of  a  single  circular  inclusion  is  modeled  by  level-1 
and  level-2  elements  under  applied  uniaxial  tension  loading.  In  the  first  case  the  En  macroscopic 
element  is  coupled  with  the  microscopically  periodic  RVE  by  asymptotic  homogenization,  while  in 
the  second  case,  the  tension  load  is  directly  applied  on  one  edge  of  the  Voronoi  cell  element  E12  model 
of  the  RVE.  The  loading  is  continued  beyond  the  level  of  inclusion  cracking  .This  is  represented 
by  the  loss  of  stress  carrying  capacity  in  the  averaged  stress-stain  plot  of  figure  4.14a.  Though  the 
response  is  similar  in  the  initial  elastic  phase,  it  diverges  with  the  onset  of  plastic  flow.  The  damage 
occurs  earlier  in  the  level-1  element  due  to  the  additional  constraint  of  periodicity.  The  contour 
plots  of  inclusion  maximum  principal  stress,  normalized  by  the  critical  stress  are  compared 
at  a  macroscopic  strain  2.1%  in  figure  4.14.  This  is  just  before  cracking  by  the  level-1  model.  A 
considerably  lower  stress  level  is  seen  for  the  level-2  model.  Such  over-estimation  of  stresses  with 
periodicity  boundtiry  condition  near  free  boundaries  necessitate  the  use  of  the  proposed  multi-level 
models. 

4.10.1  Effects  of  inhomogeneity  size  (Scale  effect) 

Neglecting  scale  effects,  that  reflect  the  actual  size  of  microstructural  constituents,  is  a  charac¬ 
teristic  of  homogenization  with  assumptions  of  statistical  periodicity  of  a  vanishingly  small  RVE. 
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While  elements  Eio  and  Ea  conform  to  this  restriction,  the  level-2  elements  E12  in  this  work  depict 
the  actual  size  of  the  microstructure  through  the  multi-level  coupling  and  hence  scale  effects  prevail. 

In  this  example  the  effect  of  particle  or  void  size  on  the  evolution  of  damage  is  investigated. 
The  two  different  microstructural  RVE’s  considered  in  figure  4.15b  (i  and  ii)  have  identical  dis¬ 
tribution  (square  edge)  and  same  particle  or  void  area  fraction  of  Vf  =  20%.  But  the  RVE  sizes 
are  different  in  that,  the  size  of  the  smaller  RVE  (i)  is  0.5mm  x  0.5mm  while  that  of  the  larger 
RVE(ii)  is  1mm  x  1mm.  The  macroscopic  structure  is  a  square  plate  with  a  square  hole,  for  which 
the  initial  level-0  mesh  with  dimensions  is  shown  in  figure  4.15a.  Only  a  quarter  of  the  plate  is 
modeled  from  symmetry  considerations.  The  smallest  size  of  macroscopic  element  allowed  in  this 
analysis  by  h  -  adaptation  is  set  to  1mm.  Thus  the  E12  elements  contain  only  one  element  in  the 
VCFEM  model  for  the  leirger  RVE(ii)  but  four  elements  in  the  VCFEM  model  for  the  smaller  RVE 
(i).  The  material  properties  for  both  composite  and  porous  materials  are  as  follows. 

Aluminum  matrix  (Elastic-Plastic):  Young’s  Modulus  (£'nj)=68.3  GPa,  Poisson  Ratio  (i/rn)=0'30. 
Initial  Yield  Stress  {Yq)=  55  MPa,  and  Post  Yield  hardening  law:  Uggy  =  Yo  +  2.08e?'„ 

Alumina  particle  (Elastic):  Young’s  Modulus  {Ec)=  344.5GPa.  Poisson  Ratio  {i/c)=  0.26,  Critical 
particle  cracking  stress  (cr"' )=0.3GPa 

The  load  is  applied  in  20  equal  increments  upto  a  total  displacement  of  1  mm  as  shown  in  figure 
4.15a. 

For  the  composite  microstructures,  figure  4.16  depicts  the  evolved  macroscopic  and  level-2  com¬ 
putational  domains  at  the  end  of  loading  stages.  The  level-0,  level-1,  level-2  elements  are  shown  in 
white,  grey  and  black  respectively  for  the  adapted  macroscopic  mesh  in  4.16a  and  b.  The  evolution 
of  the  levels  and  mesh  with  h-adaptation  is  provided  in  table  3.  The  effect  of  the  microstructure  size 
becomes  more  pronounced  with  increasing  deformation.  A  larger  level-2  domain  (29  macroscopic 
elements)  with  a  smaller  level-0  domain  is  evidenced  for  the  smaller  RVE  (i).  The  effect  of  RVE  size 
on  the  pattern  of  particle  cracking  is  very  significant.  The  level-2  region  shows  24  cracked  particles 
for  the  RVE  (i)  as  opposed  to  6  cracked  particles  for  the  RVE(ii).  The  path  of  evolution  of  cracked 
particles  is  quite  different  for  the  two  models.  For  the  RVE(ii),  the  aggregate  microscopic  cracks 
is  found  to  propagate  perpendicular  to  the  macroscopic  loading  direction  and  all  the  microcracks 
have  the  same  orientation  as  the  chain  or  macrocrack.  For  the  smaller  ptirticles  in  RVE  (i),  the 
chain  of  microcracks  or  the  macrocrack  is  observed  to  move  at  45°  to  the  loading  direction  with 
individual  microcracks  predominantly  at  90°  to  the  loading  direction. 

The  contour  plots  of  mticroscopic  and  microscopic  plastic  strain  distribution  at  the  final  loading 
stage  are  shown  in  figures  4.17  and  4.18  for  the  two  microstructures.  The  macroscopic  strain  distri¬ 
bution  for  both  models  shows  higher  strain  concentration  at  the  corner  of  square  hole  with  increased 
loading.  While  the  macroscopic  strains  aie  not  very  different  for  the  two  RVE’s,  significantly  larger 
local  plastic  straiins  are  seen  in  the  level-2  microstructure  of  RVE  (i).  A  better  jcepresentation  of 
this  difference  is  seen  in  the  graphs  of  figures  4.19  and  4.20.  The  macroscopic  (averaged)  stress  En 
history  at  the  critical  corner  in  figure  4.19  does  not  indicate  significant  difference  and  hence  exhibits 
little  scale  effect.  The  stress  drops  in  this  figure  correspond  to  particle  damage.  The  histogram 
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of  the  evolution  of  particle  cracking  however  shows  a  considerably  different  pattern  and  a  much 
higher  rate  of  cracking  for  RVE  (i). 

The  same  problem  is  solved  for  porous  microstructures  with  the  two  sized  RVE's,  but  with  a 
total  displacement  of  0.4mm.  The  evolved  levels  and  meshes  in  the  computational  models  at  the 
end  of  loading  are  shown  in  figures  4.21a  and  b  and  the  corresponding  level-2  microstructures  in 
figures  4.21c  and  d.  The  evolution  of  the  computational  domain  is  also  charted  in  table  4.  It  is 
interesting  to  note  that  the  difference  in  response  for  the  two  RVE’s  is  insignificant  this  case.  This 
may  be  attributed  to  lack  of  matrix  damage  or  softening  in  the  model.  The  plastically  hardening 
matrix  does  not  trigger  adequate  difference  in  the  adaptation  criteria  as  the  particle  cracking  does 
for  composites.  A  larger  level-1  domain  opens  up  with  the  adaptation  criteria  for  elements  which 
appear  to  be  on  the  verge  of  strain  localization,  but  subsequently  do  not  make  the  transition  to 
level-2.  The  contour  plots  of  macroscopic  and  microscopic  plastic  strain  distribution  in  figures 
4.22  and  4.23,  show  similar  macroscopic  strain  distributions,  but  higher  strains  for  the  microscopic 
model  with  smaller  porosity  RVE  (i),  due  to.  the  proximity  of  voids.  This  again  shows  the  scale 
effect  on  the  solution. 

4.10.2  Homogenization  vs.  multi-level  simulation 

The  effect  introducing  level-2  elements  on  both  macroscopic  and  microscopic  response  is  studied 
by  compeu-ing  a  pure  level- 1  simulation  of  the  square  composite  plate  in  the  previous  example  with 
a  multi-level  simulation.  The  results  shown  are  for  the  larger  particles  in  RVE(ii).  The  figure  ?? 
shows  the  microstructure  necir  the  inside  corner  by  the  two  models  at  the  end  of  loading.  The 
boxed  RVE's  in  figure  ??a  symbolize  their  periodic  repetition.  The  periodicity  constraint  results  in 
a  considerably  large  portion  of  the  microstructure  being  damaged  for  the  homogenized  simulation. 
The  direction  of  the  damage  evolution  indicated  by  the  homogenized  model  is  also  different  from 
the  level-2  simulations.  The  stress  along  the  section  A-B  is  plotted  in  figure  4.25  to  evaluate 
the  effect  of  homogenization  on  stress  concentrations  near  the  corner  and  free  edge.  The  two 
models  behave  similarly  upto  the  neighborhood  of  the  corner.  While  the  multi-level  model  predicts 
a  higher  peak  near  the  corner,  it  subsides  considerably  to  meet  the  traction  free  edge  conditions. 
The  corresponding  microscopic  level-2  stress  variations  for  the  multi-level  model  are  shown  in  the 
inset. 


4.10.3  Effect  of  heterogeneity  distribution  and  shape 

To  illustrate  the  influence  of  particle  distribution  on  the  macro-microscopic  damage  response,  two 
RVE’s  are  selected  with  same  volume  fraction  (20%),  .size  (1.0mm)  and  number  of  particles  (25). 
One  has  a  hardcore  distribution  (figure  4.15b  (iii)).  which  is  a  random  distribution  with  a  minimum 
permissible  distance  between  pauticles,  while  the  other  lias  one  cluster  in  a  hardcore  background. 
Proximity  of  particles  within  the  cluster  is  known  to  initiate  damage  in  discrete "microstructures. 
The  starting  macroscopic  mesh  is  the  same  as  in  the  previous  examples  and  a  total  displacement 
of  0.55mm  is  applied  on  the  edge  in  equal  increments.  The  smallest  allowed  size  of  macroscopic 
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elements  by  h-2idaptation  is  set  to  1mm  such  that  each  E12  element  consists  of  one  RVE. 

The  evolved  macroscopic  models  for  the  two  RVE’s,  at  the  finish  of  the  loading  cycle,  are  shown 
in  figures  4.26a  and  b.  The  levei-2  region  (black)  for  the  clustered  RVE  is  larger  than  that  for 
the  hardcore  RVE.  Within  the  E12  elements,  only  one  element  for  the  hardcore  distribution  expe¬ 
riences  particle  damage  as  shown  in  figure  4.26d.  However,  several  E12  elements  for  the  clustered 
microstructure  exhibit  particle  damage,  mainly  within  the  cluster  (figure  4.26e).  While  the  macro¬ 
scopic  averaged  plastic  strains  show  very  little  difference  for  the  different  microstructures  in  figures 
4.27a  and  4.28a,  the  microscopic  plots  in  figures  4.27b  and  4.28b  clearly  depict  the  influence  of 
distribution.  Much  higher  levels  of  effective  plastic  strain  values  are  observed  within  the  cluster, 
compared  to  significantly  lower  levels  in  the  hardcore  RVE.  Figure  4.30  shows  the  and  evolution 
of  macroscopic  stress  Eix  at  ^the  corner  of  the  square  hole  and  the  number  of  damaged  particles 
as  a  function  of  straining.  The  stress  drops  to  lower  values  for  the  clustered  RVE  due  to  a  larger 
damage  microstructure.  More  than  twice  the  number  of  particles  are  damaged  for  the  clustered 
case  as  shown  in  the  histogram. 

To  investigate  the  influence  of  shape,  a  RVE  (see  figure  4.15b  (iv))  with  the  same  volume  fraction 
(20%),  size  (1.0mm)  and  number  (25)  is  considered.  The  particles  are  elliptical  with  aspect  ratio 
3.0  and  randomly  distributed  and  oriented.  The  evolved  macroscopic  model  in  figures  4.26c  shows  a 
larger  level-2  region  compared  with  other  two  microstructures,  with  several  Ei2  elements  exhibiting 
particle  damage.  The  much  larger  number  of  cracked  particles  is  also  observed  from  the  histogram 
of  figure  4.30b.  For  this  case,  both  macroscopic  and  microscopic  plastic  strains  in  figure  4.29  are 
also  considerably  larger.  The  macroscopic  stress  En  shows  a  larger  drop  due  to  the  increased 
damage  in  the  microstructure. 

4.10.4  A  heterogeneous  plate  with  a  macroscopic  holes 

A  different  macroscopic  domain,  viz.  a  plate  with  periodically  repeated  square  diagonal  array  of 
circular  holes  is  considered  in  this  final  example.  The  plate  is  incrementally  loaded  using  prescribed 
displacement  on  the  top  and  bottom  surfaces  to  a  total  extension  of  0.15  mm.  Due  to  periodicity 
and  symmetry,  only  a  part  of  domain  is  modeled  as  shown  in  figure  4.31.  The  radius  of  the  circular 
holes  are  50mm  for  the  100  mm  x  100  mm  square  plate  as  shown  in  figure  4.31.  The  microstructural 
RVE  is  a  20%  volume  fraction  0.4  mm  x  0.4  mm  square  region  with  a  single  circular  particle.  The 
same  material  properties  as  in  the  previous  example  are  used  with  the  only  exception  being  that 
the  critical  pzu-ticle  cracking  stress  =  0.2  GPa. 

The  adapted  multilevel  computational  domain  is  shown  in  figure  4.32  and  the  number  of  el¬ 
ements  in  each  levels  with  increasing  are  tabulated  in  table  5.  The  level-2  elements  are  created 
along  a  clear  path  connecting  the  holes  due  to  localization  by  particle  cracking.  Extended  portions 
of  the  damaged  microstructure  in  level-2  regions  are  shown  in  figure  4.32.  The  macroscopic  plastic 
strain  contour  in  figure  4.33a  gives  an  indication  of  ‘hot  spots’  of  evolving  damage  near  the  central 
region.  The  microscopic  strain  plots  in  figure  4.33b  show  a  large  fraction  of  particles  cracked  and 
may  be  interpreted  as  the  initiation  of  localization  to  cause  failure  between  the  holes. 
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4.11  An  Example  on  Convergence  of  the  Multi-level  Method 

To  provide  convergence  chairacteristics  of  the  multi-level  model,  an  example  problem  of  an  elastic 
fibrous  composite  with  a  free-edge  and  loaded  in  the  out  of  plane  direction,  is  considered.  This 
classical  problem  was  proposed  by  Pagano  and  Rybicki  (1974)  to  demonstrate  the  limitations  of 
effective  modulus  or  homogenization  theory  in  predicting  stress  states  in  laminated  composites, 
especially  near  the  free-edge.  The  problem  consists  of  a  composite  cross-section  as  shown  in  figure 
(4.34a).  The  upper  half  of  the  cross-section  consists  of  n  periodic  rows  of  aligned  cylindrical  fibers, 
arranged  in  a  square  array,  while  the  bottom  half  iff  the  homogeneous  matrix  material.  The  ratio 
of  fiber  radius  to  edge  length  in  the  local  RVE  of  figure  (4.34b)  is  j  =  0.3,  corresponding  to  a 
local  volume  fraction  of  28.2%.  The  body  is  subjected  to  a  generalized  plane  strain  loading  with 
out-of-plane  loading,  prescribed  as  =  1.  Due  to  symmetry  in  the  xz  and  yz  planes  only  one 
quarter  of  the  laminate  is  modeled  (figure  4.34a).  Symmetric  boundary  conditions  are  employed 
on  the  surfaces  x  =  0  and  y  =  0,  and  the  top  and  right  surfaces  axe  assumed  to  be  traction  free. 
The  material  properties  for  the  boron  fiber  and  epoxy  matrbc  are  prescribed  as:  Eboron  =  60x10® 
psi,  I'boTon  —  0-2,  Eepoxy  —  0.5x10®  psi  and  Uepoxy  =  0.34. 

4.11.1  Microscale  VCFEM  Simulations 

To  establish  convergence  of  solutions  at  the  micro-scale,  the  adaptive  Voronoi  cell  finite  element 
model  (AVCFEM)  developed  in  Moorthy  and  Ghosh  (1999)  is  coupled  with  the  multi-level  com¬ 
putational  model.  In  AVCFEM,  two  error  measures  are  introduced  to  measure  the  quality  of  the 
solution.  They  are  (i)  the  traction  reciprocity  error,  derived  ai-posteriori  from  solved  traction 
discontinuity  at  the  element  boundary  and  matrix- inclusion/ void  interface,  and  (ii)  the  error  in 
kinematic  relationships  that  is  equated  to  an  error  in  the  strain  energy,  in  each  of  the  element  con¬ 
stituent  phases.  Displacement  adaptations  that  minimize  traction  discontinuity,  are  implemented 
to  approach  ‘optimal’  directions  for  displacement  fields  on  the  boundary /interface.  These  direc¬ 
tions  optimize  the  virtual  work  due  to  traction  discontinuity  and  are  obtained  as  components  of 
the  traction  discontinuity  in  directions,  orthogonal  to  the  original  displacement  field.  It  is  achieved 
through  h-adaptation  by  adding  displacement  degrees  of  freedom  in  the  x-  or  y-  directions  on  the 
boundary/interface.  This  is  followed  by  polynomial  or  spectral  p-enrichment  of  boundary  element 
interpolation  functions  [L].  The  error  in  kinematic  relations  may  be  minimized  by  enhancing  the 
stress  functions  by  polynomial  enrichment  or  ®"’’p-adaptation  in  the  matrix  and  inclusion  phases. 
Since  the  enrichments  are  carried  out  simultaneously  for  all  .elements,  the  adaptation  has  elements 
of  minimizing  both  local  and  pollution  errors  in  the  microstructure.  Numerical  analyses  of  several 
microstructures  with  different  distributions,  sizes  and  shapes  of  heterogeneities  have  shown  a  supe¬ 
rior  convergence  rate  with  these  adaptations. 

The  adaptive  VCFEM  of  Moorthy  and  Ghosh  (1999)  is  incorporated  to  reduce  error  in  mi- 
crostructural  simulations,  especiadly  in  level-2  elements.  To  illustrate  the  effect  of  &icroscale  adap¬ 
tation,  a  simplified  form  of  the  problem  posed  above,  with  1  row  of  4  fibers  is  considered.  The 
mesh  consists  of  4  Voronoi  cell  elements  as  shown  in  figure  4.35a.  The  matrix  stresses  in  the  hybrid 
formulation  are  constructed  from  am  Airy’s  stress  function  consisting  of  a  12  term  4-th  order  poly- 


75 


nomial  expansion  and  a  36  term  reciprocal  function  (see  Moorthy  aind  Ghosh  (1996,1999)).  The 
inclusion  stresses  are  constructed  using  a  25  term,  6-th  order  polynomial  function.  Displacement 
fields  corresponding  to  equation  (4.22)  are  constructed  with  linem  shape  functions  for  element 
boundaries  and  quadratic  shape  functions  for  curved  interface  elements.  The  sequence  of  euiapta- 
tions  consists  of  two  cycles  of  h-  followed  by  a  cycle  of  p-  adaptation  of  displacement  degrees  of 
freedom,  followed  by  a  cycle  of ''"’’p-adaptation  of  the  stress  functions.  The  adapted  mesh,  showing 
added  displacement  degrees  of  freedom,  is  illustrated  in  figure  (4.35a).  The  pre-adaptation  nodes 
are  marked  with  a  •,  the  x  direction  nodal  adaptations  are  marked  with  a  □  while  those  in  the  y  _  . 
direction  are  shown  with  A.  The  convergence  of  the  VCFEM  solutions  is  shown  in  figure  (4.35b),  " 
only  for  the  average  traction  reciprocity  error  or  {A.T.R.E.)  as  a  function  jo£\ he  total  degrees  of., 
freedom.  The  A.T.R.E.  is  calculated  as  (see  Moorthy  and  Ghosh  (1999)): 


A.T.R.E. 


fNs  ~E  pi' 

-.eE  =  \  ^  Ae,  =  l  gr 


(4.48) 


where  Ne^  correspond  to  the  total  number  of  segments  and  the  traction  reciprocity 

error  on  the  element  boundciry  and  interface  elements  respectively.  The  degrees  of  freedom  (D.O.F.) 
correspond  to  the  sum  of  the  total  number  of  nodal  degrees  of  freedom  and  the  number  of  0  stress 
parameters,  i.e.  D.O.F.  =  2*  N^odes  +  ^0-  Table  6  and  figure  4.35b  provide  numerical  details  of 
A.T.R.E  with  added  degrees  of  freedom  due  to  each  adaptation,  on  the  boundary  and  interface 
only.  Discrete  points  in  the  figme  correspond  to  the  different  stages  of  adaptation.  The  first  two 
drops  are  for  the  two  consecutive  cycles  of  /i-adaptation  and  the  third  for  p-zwiaptation.  The  trac¬ 
tion  reciprocity  error  reduces  rather  drastically  in  the  /i-adaptation  cycles,  i.e.  a  90%  change  in 
A.T.R.E.  is  obtained  with  a  61%  increase  in  D.O.F..  With  the  subsequent  p-  adaptation,  and 
additional  4%  change  in  A.T.R.E.  is  obtained  with  a  further  13%  increase  in  D.O.F..  Since  very 
little  A.T.R.E.  reduction  is  gained  beyond  two  cycles  of  boundeiry  /i-adaptations,  only  these  are 
employed  for  the  multi-scale  simulations  of  this  section. 


A  comparison  study  is  made  with  the  solutions  in  Pagano  and  Rybicki  (1974)  (the  same  results 
arte  also  generated  by  the  commercial  code  ANSYS).  The  ANSYS  mesh,  ^^r' which  convergence 
is  thieved  has  4230  QUAD4  qleinents  and  4352  nodes.  The  microscopic  jn:piane  stress  cTyy  plots 
along  horizontal  sections  y  =  h  are  compared  for  the  unadapted  VCFEM,  adapted  VCFEM  and  the 
ANSYS  model  in  figure  4.36.  It  should  be  noted  that  the  transverse  stress  ayy  is  approximately  two 
orders  lower  compared  to  the  leading  order  stress  and  hence  convergence  is  harder  to  achieve. 
Results  of  the  adapted  VCFEM  results  agree  very  well  with  those  of  the  converged  ANSYS  model. 
This  is  therefore  a  strong  attestation  of  the  accuracy  of  the  adaptive  VCFEM  solutions.  The  figure 
also  shows  the  singular  natmre  of  the  solution  in  the  free  edge  region  when  effective  modulus  theory 
by  homogenization  is  used.  Since  this  is  physically  unattainable,  detailed  micromechanics  solutions 
are  needed  in  this  region.  c. 
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4.11.2  Multi-level  Simulations 

To  test  convergence  of  the  overall  model  qualitatively,  the  composite  problem  with  40  (n  =  40) 
rows  of  6400  fibers  is  simulated  by  the  multi-level  code.  The  level-0/ 1  h-adaptation  criterion  is 
based  on  a  traction  jump  criterion  across  adjacent  element  boundaries,  stated  as: 


Refine  element  if  E,  >  Cl  *  Eavg  where  Eavg 


E  f^y  {m\]'^  +  [\Ty\?f^ddY 
NE  '  ‘  /ay  day 


(4.49) 

where  NE  is  the  total  number  of  level-0/ 1  elemsnts.  This  criterion  is  intended  for  signaling  and 
zooming  in  on  regions  high  stress  gradients.  For  level-0  to  level- 1,  the  transition  is  made  based  on 
element  jump  in  the  significant  stress  Cyy,  in  addition  to  the  traction  jump  condition  of  equation 
(4.49),  stated  as: 


Level  —  0  — ^  Level  —  1  if  Ej  ^  Cl* Eavg  or  ^E^  ^  Cl *^ Eavg 


where  ^Ej  = 


hvMddY 

layddY 


(4.50) 


where  C1,C2  are  prescribed  constants.  Finally,  the  level-1  to  level-2  transition  is  made  from 
observations  made  in  the  microstructural  RVE’s  of  level-1  elements.  It  is  based  on  the  ratio  of 
local  energy  density  (at  each  integration  point  of  VC  element)  to  the  average  energy  density  for 
the  entire  RVE  (aggregate  of  VC  elements).  The  procedure  followed  for  determining  this  is: 

•  Solve  the  RVE  problems  with  periodicity  boundary  conditions  eind  four  sets  of  applied  macro¬ 
scopic  strains,  viz.  (i)  en  =  1,  eyy  =  0,  exy  =  0,  ezz  =  0,  (ii)  ej*  =  0,  eyy  =  1,  exy  =  0,  ezz  =  0, 
(iii)  eix  =  0,  eyy  =  0,  exy  =  1,  Czz  =  0  and  (iv)  Cn  =  0,  eyy  =  0,  Cxy  =  0,  Czz  =  1- 

•  For  all  cases  (i-iv)  evaluate  the  energy  densities  at  integration  point  of  the  matrix  and  inclusion 

phases  as  where  Siju  is  the  compliance  tensor. 

•  Evaluate  the  energy  density  concentration  factors: 

,  _  Maximum  U^)  ^  Maximum 

“  RVE  Averaged  C'"“‘(=  C^S?*)  RVE  Averaged  t/*"^(= 


•  In  the  multi-level  problem  with  actual  strain  components,  the  transition  is  made  according 
to: 


Level  —  1  Level  —  2  if 


actual  rrmat 


Tj-mai 
^max  — 


^mat  actual  rrmat 


or 


actual  Trincl 


Uinci  > 
^max  - 


R 


incl  actual  jjind 

^  ^  nvpT 


(4.51) 


at  more  than  1%  of  all  integration  points. 


The  initial  mesh  consists  of  200  level-0  QUAD4  elements.  To  examine  the  dependence  of  results  on 
adaptation  parameters,  two  sets  of  parjimeters  viz.  (Cl  =  1.5,  C2  =  2)  and  (Cl  =  1.25,  C2  =  1.5) 
axe  experimented  with.  Figures  4.37a  and  b  show  the  h-adapted  level-O/level-l/level-2  meshes  with 
the  diflferent  adaptation  parameters.  The  first  set  leads  to  391  level-0  elements,  0  level-1  elements. 
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6  transition  elements  and  4  level-2  elements,  while  the  second  set  has  475  level-0  elements,  4  level- 
1  elements,  9  transition  elements  and  6  level-2  elements.  The  microstructural  RVE  or  unit  cell 
is  assumed  to  be  of  dimension  I  =  ^  and  the  level-2  elements  are  assumed  to  contain  a  single 
unit  cell.  The  same  problem  is  also  solved  using  the  commercial  code  ANSYS  with  a  mesh  of 
30.000  elements,  for  which  a  3  x  3  array  of  fibers  near  the  free-edge  and  laminate  interface  are 
modeled  explicitly  (see  4.37c).  In  figure  4.38a  the  macroscopic  response  of  level-0  elements,  with 
homogenized  effective  moduli,  is  plotted.  The  figure  shows  the  effect  of  mesh  dependence  in  the 
solutions  near  the  free-edge.  With  increasing  cycles  of  adaptations  the  artificial  stress  singularity 
moves  towards  the  free  edge  and  gains  intensity.  In  figure  4.38b  the  mi-'^oscopic  level-2  stresses  are 
plotted  in  the  vicinity  of  the  free  edge  (indicated  as  L2  in  figure  4.38a),  It  is  clearly  seen  that  are 
artificial  singularity  of  homogenized  solutions  no  longer  exist  and  there  is  no  mesh  dependence  near 
the  free  edge.  The  two  multi-level  meshes  of  figures  4.38  a  and  b  and  the  ANSYS  model  all  produce 
similar  results.  Thus,  through  this  example,  the  convergence  characteristics  of  the  multi-level  model 
is  qualitatively  verified  and  confidence  is  gained  about  its  application. 

4.12  Discussions  and  Conclusion 

Adaptivity  in  the  computational  modules  for  multiple  scale  problems  entails  minimizing  two  types 
of  errors,  viz.  the  discretization  error  and  the  modeling  error.  In  this  work  an  adaptive  multi-level 
method  is  proposed  to  primarily  focus  on  reducing  the  modeling  error  and  predicting  the  evolution 
of  stresses,  strains  and  damage  at  the  structural  and  microstructural  scales.  The  microstructural 
analysis  is  conducted  with  the  Voronoi  cell  finite  element  model  (VCFEM)  for  elastic-plastic  con¬ 
stituents  with  particle  cracking.  VCFEM  allows  for  continuous  change  in  element  topology  due  to 
progressive  damage  with  high  accuracy  as  shown  in  Moorthy  and  Ghosh  (1999)  and  Ghosh  and 
Moorthy  (1998).  The  efficiency  of  the  method,  due  to  embedding  micromechanics  in  FEM  formula¬ 
tion,  makes  modeling  of  large  microstructural  regions  relatively  easy.  A  conventional  displacement 
based  elastic-plastic  FEM  code  is  developed  for  macroscopic  analysis.  Adaptive  mesh  refinement 
and  level  transition  strategies  are  developed  to  create  a  hierarchy  of  computational  sub-domains 
with  varying  resolutions.  This  differentiates  between  non-critical  and  critical  regions  and  helps  in 

increasing  the  efficiency  of  computations  by  preferential  "zoom  in’. 

Three  levels  of  hierarchy,  viz.  level-0,  level-1  and  level-2,  evolve  in  the  multi-scale  model  with 
progressive  deformation.  A  piecewise  continuous  elastic-plastic  constitutive  law  is  developed  for 
level-0  simulations.  The  specialty  of  this  model  is  that  since  it  is  developed  from  rigorous  mi¬ 
cromechanical  simulations  with  precise  material  morphology,  it  is  sensitive  to  variations  in  the 
microstructural  distribution.  The  constitutive  relation  leads  to  very  high  efficiency  in  simulations 
when  compared  with  two-scale  analysis  by  homogenization  (e.g.  Guedes  and  Kikuchi  (1991)). 

Level-0  simulations  are  accompanied  by  mesh-refinement  using  h-adaptation  techniques,  to  re¬ 
duce  discretization  error  in  the  computational  model  and  zoom  in  on  regions  of  evolving  localization 
due  to  microscopic  non-homogeneity.  The  criteria  for  h-adaptation  are  based  on  gradients  or  jumps 
in  plastic  work  or  stresses.  While  these  error  criteria  are  effective  in  equi-distributing  discretization 
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errors  where  solutions  are  relatively  regular,  they  are  also  helpful  in  identifying  regions  of  localiza¬ 
tion  or  singularity  due  to  steep  gradients  in  the  solution  variables.  It  should  however  be  mentioned 
that  the  proposed  refinement  schemes  and  error  criteria  are  by  no  means  optimal.  For  example, 
it  has  been  shown  in  Guo  and  Babuska  (1986),  Rank  and  Babuska  (1987)  the  h  —  p  adaptation, 
which  combines  local  mesh  refinement  with  increase  of  the  polynomial  order  of  interpolants,  enjoy 
exponential  rate  of  convergence  and  excellent  accuracy  especially  for  problems  with  singularity. 
Additionally,  it  has  been  demonstrated  in  Oden  and  Feng  (1996)  and  Babuska  et.  al.  (1994), 
that  local  error  estimation  techniques  are  incapable  of  detecting  pollution  error  particularly  when 
singularities  are  present.  A  robust  mesh  adaptation  procedure  for  reducing  errors  should  therefore 
incorporate  pollution  error  reduction  criteria,  £is  well  as  h  —  p  adaptation.  Much  of  the  ongoing 
and  future  efforts  by  the  authors  are  focussed  on  introducing  these  techniques. 

When  imminent  damage  and  localization  are  sensed  by  the  code,  the  level-0  elements  automat¬ 
ically  switch  to  level- 1  elements,  which  use  computations  both  at  the  macroscopic  and  microscopic 
scales.  The  criteria  for  signaling  such  transition  are  nonunique  and  a  number  of  them  are  con¬ 
sidered  in  this  study.  Among  those  that  yielded  effective  results  are;  (i)  when  the  homogenized 
constitutive  relation  has  reached  its  limit;  (ii)  when  the  maximum  macroscopic  principal  stress  or 
hydrostatic  stress  exceeds  a  certain  fraction  of  the  microscopic  fracture  stress  or  (iii)  when  the 
macroscopic  effective  plastic  strain  or  the  dilatation  exceeds  certain  pre-determined  values.  With 
the  rise  in  local  gradients  of  macroscopic  veiriables  or  with  microstructural  damage,  the  pre-assumed 
representative  volume  element  in  the  microstructure  are  no  longer  effective  and  a  shift  to  complete 
microscopic  simulations  is  necessary.  Extended  portions  of  the  microstructure  are  directly, mod¬ 
eled  by  VCFEM  in  these  level-2  elements.  An  adaptive  VCFE  model  (AVCFEM)  is  used  to  study 
convergence  of  the  microstructural  analysis  model.  In  AVCFEM,  traction  reciprocity  error  and 
error  in  kinematic  relationships  are  optimized  through  boundary/interface  enrichment  in  the  form 
of  additional  displacement  degrees  of  freedom  and  interpolation  orders,  as  well  as  through  stress 
function  enhancement.  In  Moorthy  and  Ghosh  (1999),  excellent  convergence  rate  of  AVCFEM  has 
been  established  for  elastic  and  elasto-plastic  microstructural  analysis.  For  the  problem  considered 
in  this  work,  two  cycles  of  boundary  h-  adaptation  is  found  to  achieve  desirable  results  in  the 
level-2  elements. 

Several  numerical  examples  are  conducted  with  the  multi-level  model  to  examine  the  efface  of 
various  microstructural  morphologies  on  the  multi-scale  response  of  composite  and  porous  struc¬ 
tural  components.  Specifically  scale  effects,  effects  of  microstructural  distribution  and  shapes  and 
structural  geometry,  on  the  mechanical  and  damage  response  are  investigated.  In  conclusions, 
the  multi-level  model  has  addressed  a  number  of  difficult  issues  in  solving  multi-scale  problems. 
However  a  number  of  avenues  for  further  model  enhancements  remain.  For  example,  the  Newton- 
Raphson  iteration  scheme  used  for  nonlinear  problems  struggle  for  convergence,  as  more  and  more 
level-2  elements  with  damage  emerge.  The  migration  of  boundaries  between  different  levels  within 
each  iteration  loop  has  an  eidverse  effect  on  convergence  as  well.  Also  for  the  levfftl-l  elements,  an 
iterative  loop  is  necessary  for  convergence  of  steps  (a),  (b)  and  (c)  in  section  4.6.  In  addition, 
the  issues  of  adaptivity  mentioned  above  remain.  These  issues  are  currently  being  dealt  with  for 
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maJcing  this  an  effective  tool  in  the  prediction  of  structural  failure. 
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Figure  4.1:  A  heterogeneous  structure  showing  various  scales:  (a)  The  structure  at  the  macroscopic 
scale  of  applied  loads  (b)  A  representative  volume  element  (RVE)  at  the  microscopic  scale  with 
the  VCFE  model  and  (c)  A  Voronoi  cell  element  at  the  scale  of  a  single  heterogeneity  or  basic 
structural  element. 
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Figure  4.2;  A  hierarchical  multi-level  computational  domain;  Level-0  for  macroscale  continuum 
modeling  (a)  Level-1  for  coupled  macro-microscopic  (R\  E)  modeling  with  asymptotic  homogeniza¬ 
tion  and  (b)  Level-2  for  pure  microscopic  modeling 


Figure  4.3:  Flow  chart  of  the  sequence  of  operations  in  the  multiple-scale  model. 


87 


(C,SiX,Eyj,,Eiy) 


Figure  4.4:  A  en  -  Cyy  -  e^y  strain  sub-space,  discretized  into  ci  bic  elements  for  interpolating 
constitutive  model  parameters.  The  nodal  values  of  stresses  and  plastic  work  are  numerically 
generated. 


Figure  4.5:  Variation  of  the  yield  function  coefficient  C(exi,  Cyy,  Ciy,  Wp),  (a)  as  a  function  of  0  or 
normal  strains  for  different  microstructures,  (b)  as  a  function  of  plastic  work  Wp. 
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Figure  4.7:  Comparison  of  the  stress-strain  results  in  the  coinposite  for  the  uniaxial  loading  test 
by  (i)  macroscopic  analysis  with  the  homogenized  constitutive  model  amd  (ii)  two-scale  analysis 
with  asymptotic  homogenization  model  for:  (a)  Cl.  (b)  C2,  (c)  C3,  (d)  C4,  (e)  C5  and  (f)  C6 
microstructures. 
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Figure  4.10:  Contour  plot  of  the  microscopic  plastic  .strain  in  the  voided  microstructures  under 
pure  shear  loading  condition,  for  (a)  VI  (b)  V2  microstructures. 
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Figure  4.11:  Comparison  of  the  stress-strain  results  in  the  porous  material  for  (a)  bi-2ixial  ten¬ 
sion  for  VI  microstructure  (b)  bi-axial  tension  for  V3  niicrostructure,  (c)  uniaxial  tension  for  VI 
microstructure  and  (d)  uniaxial  tension  for  V3  microstructure. 
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Figure  4.7:  Compzurison  of  the  stress-strain  results  in  the  composite  for  the  uniaxial  loading  test 
by  (i)  macroscopic  analysis  with  the  homogenized  constitutive  model  and  (ii)  two-scale  analysis 
with  asymptotic  homogenization  model  for:  (a)  Cl,  (b)  C2,  (c)  C3,  (d)  C4,  (e)  C5  and  (f)  C6 
microstructures. 
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Figure  4.10;  Contour  plot  of  the  microscopic  plastic  strain  in  the  voided  microstructures  under 
pure  shear  loading  condition,  for  (a)  VI  (b)  V2  microstructures. 
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Figure  4.11:  Comparison  of  the  stress-strain  results  in  the  porous  material  for  (a)  bi-axial  ten¬ 
sion  for  VI  microstructure  (b)  bi-axial  tension  for  V3  inicrostructure,  (c)  uniaxial  tension  for  VI 
microstructure  and  (d)  uniaxial  tension  for  V3  microstructure. 
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Figure  4.12;  Creating  extended  microstructure;  (a)  a  mesh  of  macroscopic  elements  with  an  under¬ 
lying  microstructure  of  repeated  RVEs,  (b)  the  extended  microstructure  by  Voronoi  tessellation. 
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Figure  4.13:  Interfaces  between  the  levelO/1  Eio  and  En  elements,  transition  Etr  elements  and  the 
level-2  E(2  elements. 
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Figure  4.14:  (a)  Macroscopic  axial  stress-strain  plots  for  a  single  RVE  in  tension  modeled  by  level- 1 
and  level-2  elements;  Maximum  principal  stress  distribution  in  the  inclusion  by  (b)  level-2  element 
and  (c)  level-1  element.  ^ 
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Figure  4.15:  (a)  A  quarter  of  the  level-0  starting  macroscopic  model  of  a  square  plate  with  a  square 
hole,  (b)  four  different  microstructural  RVE’s  for  the  [)late  model.  ^ 
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Figure  4.16;  Macroscopic  three  level  evolved  mesh  for  the  composite  material  at  the  end  of  the 
loading  cycle  for  (a)  the  larger  RVE  (ii)  and  (b)  smaller  RVE  (i);  the  corresponding  level-2  region 
with  the  (c)  the  larger  RVE  (ii)  and  (d)  smaller  RVE  (i).  Level-0  is  with  white  elements,  level-1  is 
with  grey  elements  and  level-2  is  with  black  eleipgnts. 
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Figure  4.17:  Contour  plot  of  plastic  strain  ^  for  the  composite  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  larger  RVE  (ii) 
model. 
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Figure  4.18:  Contour  plot  of  plastic  strain  eP  for  the  composite  square  plate:  (a)  ^e  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  smaller  RVE  (i) 
model. 
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Figure  4.19:  The  evolution  of  Sn  at  the  corner  node  of  square  hole  for  the  two  microstructures. 
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Figure  4.20:  Histogram  of  the  number  of  damaged  particles  as  a  function  of  straining. 


Figure  4.21:  Macroscopic  three  level  evolved  mesh  for  the  porous  material  at  the  end  of  the  loading 
cycle,  for  (a)  the  larger  RVE  (ii)  and  (b)  smaller  RVE  (i);  The  corresponding  level-2  region  with 
the  (c)  the  larger  RVE  (ii)  and  (d)  smaller  RVE  (i).  Level-0  is  with  white  elements,  level-1  is  with 
grey  elements  and  level-2  is  with  black  elementq^Ol 
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Figure  4.22:  Contour  plot  of  plastic  strain  ^  for  the  porous  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  smaller  RVE  (ii) 
model. 
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Figure  4.23:  Contour  plot  of  plastic  strain  for  the  [xirous  square  plate:  (a)  ^he  macroscopic 
averaged  strain  zind  (b)  the  level-2  microscopic  strain  at  tlu'  critical  region  for  the  smaller  RVE  (i) 
model. 


102 


□ 

' 

□ 

□ 

□ 

□ 

□ 

□ 

o  o 

• 

• 

0 

o 

1 - 

• 

• 

o  o 

• 

o 

• 

• 

• 

o  o 

• 

• 

• 

• 

»  o 

• 

(a)  (b) 


Figure  4.24:  A  comparison  of  the  (ji:^y:ipstructural  evolution  near  the  inside  corner  by  (a)  'complete 


Figure  4.25:  Comparison  of  macroscopic  stress  Exi  along  section  A-B  by  level-1  and  multi-level 
analysis. 
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Figure  4.26:  Macroscopic  three  level  evolved  mesh  for  the  composite  material  at  the  end  of  the 
loading  cycle  for  (a)  hardcore  distribution  with  circular  particles,  (b)  clustered  distribution  with 
circular  particles,  and  (c)  hardcore  distribution  with  elliptical  particles;  the  corresponding  level-2 
microstructures  for  (d)  hardcore  (e)  clustered  and  (f)  elliptical  RVE’s. 
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Figure  4.27:  Contour  plot  of  plastic  strain  for  the  composite  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  heirdcore  RVE 
model. 


Figure  4.28:  Contour  plot  of  plastic  strain  for  the  composite  square  plate:  (a)  ^he  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  clustered  RVE 
model. 
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Figure  4.29:  Contour  plot  of  plastic  strain  for  the  composite  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  elliptical  RVE 
model. 
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Figure  4.30:  The  evolution  of  (a)  Ex*  at  the  corner  node  of  square  hole  and  (b)  the  number  of 
damaged  particles  for  the  three  microstructures. 
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Figure  4.31:  (a)  Finite  element  model  for  a  pa: 
(b)  the  VCFE  model  of  the  microstructural  cc 
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Figure  4.32:  Macroscopic  three  level  eyolved  mesh  for  the  composite  plate  with  two  holes  at  the 
end  of  loading.  The  level-2  regions  in  the  microstructure  are  shown  at  the  two  regions  A  and  B. 
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Figure  4.33;  Contour  plot  of  (a)  the  macroscopic  plastic  strain  ^  and  (b)  level-2  microscopic  plastic 
strain  f  at  the  region  C  of  the  microstructure  in  figure  4.32 
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Figure  4.34:  (a)  A  composite  cross-section  with  n  periodic  layers  of  circular  fibers,  (b)  a  represen¬ 
tative  volume  element  for  the  effective  medium  theory  calculations. 
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Figure  4.36:  Oyy  distribution  at  section  AA  (y=h)  for  the  composite  section  with  1  row  of  fiber 
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Figure  4.37:  (a,b)  H-adapted  multi-level  mesh  showing  level-0,  level- 1  and  level-2  for  two  different 
parameters  in  transition  criteria,  and  (c)  a  corresponding  ANSYS  mesh  with  detailed  raicrostructure 
modeled  near  the  free  edge. 
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Figure  4.38:  Convergence  of  stress  cXyy  at  the  section  A- A  ly  —  h):  (a)  macroscopic  stress  plots 
using  the  homogenized  moduli  and  level-0  elements,  (b)  iiu(Toscopic  level-2  stress  plots  near  the 
critical  free-edge  region. 
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Chapter  5 


Experimental-Computational 
Investigation  Of  Damage  Evolution  In 
Discontinuously  Reinforced 
Aluminum  Matrix  Composite 

5.1  Introduction 

The  corninercial  use  of  particle- reinforced  metal  matrix  composites  in  automotive,  aerospace  and 
other  engineering  systems  has  increased  in  the  last  few  decades  due  to  their  potentially  superior 
mechanical  properties,  as  well  as  their  ability  to  reduce  life-cycle  costs  through  enhanced  thermal 
stability  and  weight  reduction.  The  property  advantages  of  these  materials  are.  however,  often  di¬ 
minished  by  general  degradation  of  failure  properties  like  ductility  and  fracture  toughness.  Various 
experimental  and  numerical  studies  [1,  2,  3,  4,  5,  6,  7,  8,  9,  10]  have  been  conducted  to  understand 
the  influence  of  morphological  factors  such  as  volume  fraction,  size,  shape  and  spatial  distribution 
as  well  as  constituent  material  and  interface  properties  on  the  deformation  and  damage  behavior. 
These  studies  have  concluded  that  failure  mechanisms  are  highly  sensitive  to  local  reinforcement 
distribution,  morphology,  size,  intc-Tacial  strength  etc. 

Traditionally  unit  cell  models  [11,  12,  13,  14,  15]  based  on  the  finite  element  analysis  have 
been  used  to  predict  the  onset  and  growth  of  evolving  damage  in  composite  materials.  While  these 
models  provide  valuable  insights  into  the  microstructural  damage  processes,  simple  morphologies 
idealize  actual  microstructures  for  many  engineering  materials  that  bear  little  relationship  to  the 
actual  stereographic  features.  These  deficiencies  have  been  circumvented  in  [16,  8,  17],  where  com¬ 
putational  models  of  discontinuously  reinforced  materials  with  random  spatial  dispersion  have  been 
considered.  Richmond  and  coworkers  [18,  19]  have  investigated  the  effect  of  morphology  on  damage 
in  composite,  porous  aind  polycrystalline  materials  by  modeling  actual  geometries  obtained  from 
2D  micrographs.  Using  the  Voronoi  Cell  finite  element  model,  Ghosh  et.  al.  [9,  10]  have  examined 
the  effect  of  various  spatial  dispersions  and  particle  shape  and  size  on  the  damage  initiation  and 
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evolution  process  in  ductile  matrix  composites. 

Many  characterization  studies  with  2-D  microstructures  e.g.  [20,  21,  22,  23]  have  also  been 
conducted  to  understand  the  relation  between  microstructural  morphology  and  damage.  Experi¬ 
mental  research  [25,  26,  27,  28,  29]  has  however  pointed  to  the  necessity  of  examining  the  full  3D 
characteristics  for  understanding  the  damage  process.  These  studies  infer  that  2-D  assessment  can 
sometimes  be  misleading,  especially  in  the  presence  of  spatial  clustering.  Non-destructive  evalua¬ 
tion  methods,  e.g.  techniques  based  on  ultrasonics  e.g.  [30],  acoustic  emission  [4]  and  X-ray  based 
computer  tomography  (CT)  [31,  32]  have  emerged  as  potential  methods  for  studying  3D  damage. 
However  many  of  these  systems  are  thus  far  not  capable  of  achieving  spatial  resolutions  required  to 
accurately  capture  microscopic  particles  and  damage  in  the  particle  reinforced  MMC’s.  Buffiere  et. 
al  [33]  are  developing  a  CT  technology  to  yield  tomographic  images  with  a  higher  spatial  resolution. 

This  work  deals  with  a  combined  experimental-computational  approach  to  study  the  evolution 
of  microscopic  damage  to  cause  complete  material  failure  in  commercial  SiC  particle  reinforced 
aluminum  alloys  or  DRA’s.  Through  a  combination  of  2D  and  3D  characterization  and  analysis 
models,  it  is  intended  to  understand  what  aspects  of  microstructural  morphology  that  are  most 
critical  for  damage  nucleation  and  evolution.  Since  it  is  difficult  to  identify  the  microcrack  growth 
process  once  a  material  has  failed  completely,  an  interrupted  testing  technique  is  designed.  Sub¬ 
sequently,  sample  microstructures  in  the  severely  necked  region  are  microscopically  examined  in 
3D  using  a  serial  sectioning  method  discussed  in  [24.  25,  26].  Computer  simulated  equivalent  mi¬ 
crostructures  are  tessellated  into  meshes  of  2-D  and  3-D  Voronoi  cells.  Various  characterization 
functions  of  geometric  parameters  are  generated  and  a  sensitivity  analysis  is  conducted  to  explore 
the  influence  of  morphological  parameters  on  damage.  2D  characterization  functions  are  compared 
with  3D  to  evaluate  the  effectiveness  of  modeling  the  2D  micrographs.  Modeling  of  the  initia¬ 
tion  and  propagation  of  damage  is  conducted  with  Voronoi  Cell  Finite  Element  Method  (VCFEM) 
[9,  10,  34,  35].  Each  Voronoi  cell  element  may  consist  of  a  matrix  phase,  an  inclusion  phase  and 
a  crack  phase.  Damage  initiation  by  particle  cracking  is  assumed  to  follow  a  maximum  principal 
stress  based  Rankine  criterion.  The  VCFEM  for  particle  cracking  has  shown  a  significant  promise 
in  modeling  large  aggregates  of  heterogeneities.  While  the  appropriateness  of  3D  analyses  is  rec¬ 
ognized  for  this  study,  the  3D  VCFEM  (under  development)  does  not  currently  have  all  necessary 
features.  Due  to  enormous  computing  requirements  of  conventional  3D  FEM  models,  various  stud¬ 
ies  have  resorted  to  simplified  manifestations  of  complex  geometries  and  properties  e.g.  [7,  41,  15]. 
This  study  is  restricted  to  2D  in  the  form  of  VCFEM  analyses  of  section  micrographs.  Finally,  the 
effect  of  size  and  characteristic  lengths  of  representative  material  element  (RME)  on  the  extent  of 
damage  in  the  model  systems  is  also  investigated. 

5.2  Experiments  for  Damage  Assessment 

5.2.1  Interrupted  tests 

The  material  analyzed  in  this  work  is  a  discretely  reinforced  commercial  aluminum  tBat  is  fabricated 
by  a  powder  metallurgy  process  [36].  It  consists  of  extruded  commercial  X2080  aluminum  alloy 
with  15%  volume  fraction  SiC  particles.  The  X2080  matrix  has  a  nominal  alloy  composition  with 
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weight  percentages  of  3.8%  Cu,  1.8%  Mg  and  0.2%  Zr,  in  addition  to  low  impurity  contents  of  Fe 
and  Si.  The  precipitation  hardened  X2080  aluminum  alloy  system  is  naturally  aged  by  heat  treating 
for  4  hours  at  9Z0°F,  followed  by  cold  water  quench  and  aging  for  2  days  at  room  temperature. 

An  important  object  in  this  failure  study  is  to  obtain  adequate  microstructural  data  that  depict 
the  growth  of  damage  into  a  major  failure  path.  In  general,  it  is  difficult  to  identify  the  domi¬ 
nant  damage  mechanisms  and  also  the  microcrack  growth  process,  once  a  material  has  fractured 
completely.  Thus  an  interrupted  testing  technique  is  designed  where  the  load  and  deformation  are 
halted  in  the  material  instability  zone,  following  necking  but  prior  to  fracture.  The  tests  assume 
that  the  major  cracks  are  essentially  prominent  is  this  stage,  and  are  helpful  in  understanding  the 
linkage  mechanism  of  microcracks  or  particle  debonds  to  facilitate  growth  of  the  dominant' damage. 
To  initialize  the  testing,  estimates  of  the  necking  and  fracture  strains  are  first  obtained  by  observing 
the  behavior  of  a  tension  test  to  failure.  The  uniaxial  tension  tests  are  executed  on  a  MTS  810  ma¬ 
terial  system  with  a  HP  7044  X-Y  recorder  to  monitor  the  loads  and  strains,  and  the  critical  strains 
are  measured  with  a  MTS  632.11  strain  gauge  extensometer.  Following  the  initialization,  strain 
controlled  interrupted  tests  axe  carried  out,  in  which  the  specimens  are  loaded  to  the  instability 
region  before  the  load  is  stopped. 

Figure  5.1a  shows  a  typical  tension  specimen  for  the  naturally  aged  material.  Data  for  six 
specimens  of  this  material,  viz.  tl,  t2,  t3,  t4.  t5  and  t6  are  tabulated  in  table  1.  The  specimens 
tl,  t2,  t4  and  t6  are  obtained  from  the  outer  annulus  region  of  the  stock  material  while  t3  and 
t5  are  from  the  central  core  regions.  The  initialization  of  the  test  to  study  the  entire  material 
behavior  and  estimate  the  post-instability  region  is  done  with  specimens  tl  and  t2.  The  material 
load-displacement  curve  is  plotted  in  figure  5.1b,  from  which  the  necking  strain  is  obtained,  from 
the  peak  load  value.  For  the  specimen  tl,  the  test  is  conducted  at  a  strain  rate  is  e  =  5xl0”‘‘sec"^ 
and  the  necking  strain  and  fracture  strain  are  found  to  be  e„  =  9.15%  and  Cj  =  9.40%  respectively. 
The  short  instability  region  in  tl  prompts  a  reduced  strain  rate  e  =  3xl0“‘*sec"^  for  specimen  t2, 
for  which  €„  =  9.05%  and  Ci  =  9.20%. 

In  table  1.  e,  e„  and  e,  correspond  to  the  strain  rate,  the  necking  strain  and  the  interrupted 
strain  respectively.  The  interrupted  strain  coincides  with  the  fracture  strain  in  the  event  that 
fracture  precedes  the  load  stoppage.  This  is  indicated  with  F  or  I  in  the  table.  Load  interruption  is 
only  possible  for  the  specimens  t3  and  t6  due  to  the  extremely  short  post-instability  range  of  this 
material  in  comparison  with  the  resolution  of  the  loading  mechanism.  The  necking  Ltrains  for  the 
specimens  tl,  t2,  t4  and  t6  are  in  the  range  of  9.00%  ~  9.30%,  while  those  for  spec.inens  t3  and  t5 
are  in  the  9.80%  10.20%  range.  This  difference  is  possibly  due  to  gradients  created  by  the  heat 

treatment  at  different  locations  in  the  stock  material.  The  core  cools  slower  and  more  uniformly 
regions  near  the  surface.  This  results  in  the  more  uniform  microstructure  and  larger  necking  and 
fracture  strains  for  specimens  (t3  and  t5)  located  near  the  core  of  the  stock  material. 

5.2.2  Damage  examination  and  microscopic  analysis 

To  examine  the  dependence  of  microstructural  damage  on  the  local  morphology,eSerial  sectioning 
of  sample  coupons  extracted  from  the  load-interrupted  specimens  t3  and  t6,  is  invoked.  This 
method,  discussed  in  [24,  25,  26],  involves  gradual  removal  of  material  layers  to  obtain  a  series 
of  scanning  electron/optical  micrographs,  representing  sections  of  a  microstructure.  It  is  a  very 
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effective  method  for  reconstructing  3D  microstructures  from  a  series  of  2D  sections  of  psirticulate 
reinforced  composites,  requiring  a  resolution  of  few  microns.  Prior  to  sectioning,  locations  are 
selected  in  hgure  5.1a  for  cutting  out  the  sample  coupons.  X-rays  and  acoustic  microscopy  with 
a  AEROTECH  UNIDEX  11  acoustic  microscope  with  a  resolution  of  about  bOfim  are  used  in 
this  process  to  detect  regions  that  contain  major  crack  paths.  Polished  surfaces  of  these  extracted 
samples  are  then  examined  by  a  Nikon  optical  microscope  for  major  damage  sites.  For  the  specimen 
t3,  shorter  cracks  passing  through  2  ~  3  particles  at  most  are  found.  However,  for  the  specimen 
t6,  a  larger  crack  passing  through  5  ~  6  particles  is  identified,  and  is  consequently  chosen  for 
analysis.  Coupons  of  approximate  size  6  mm  x  6  mm  x  6  mm  are  subsequently  prepared  for 
the  serial  sectioning  operation  to  sequentially  expose  parallel  sections  of  the  microstructure.  As 
discussed  in  [24,  25,  26],  parallel  layers  in  a  direction  perpendicular  to  the  straining  direction  (see 
figure  5.1a)  are  removed  using  a  precision  dimple  grinder.  The  depth  of  material  removal  per  step 
is  selected  such  that  each  particle  is  sectioned  at  least  once,  ensuring  that  all  particles  of  interest 
are  adequately  captured  in  the  micrographs.  For  the  DRA  considered,  the  particle  size  range  is 
approximately  3-25  fxm,  with  an  average  size  of  ~  9.2  ^m  and  the  standard  deviation  is  3.891  nm. 
The  section  to  section  step  size  is  chosen  to  be  2  /xm,  corresponding  to  a  total  traversed  thickness 
of  36/jm  for  18  sections.  Two  typical  micrographs  showing  damage  are  depicted  in  figures  5.2, 
for  which  the  horizontal  corresponds  to  the  loading  direction.  The  micrographs  are  then  serially 
stacked  using  a  graphic  software  [37]  to  yield  3D  microstructures  as  shown  in  figure  5.3a.  The 
precise  3D  location,  shape,  size  and  orientation  of  each  particle  can  be  obtained  at  a  fairly  high 
resolution  by  this  method. 

5.2.3  Major  observations 

The  micrographs  of  serial  sections  3  and  5  in  figure  5.2,  perpendicular  to  the  middle  plane  of  the 
tensile  specimen,  provide  important  information  on  the  evolution  of  the  dominant  damage  path 
in  the  material.  A  dominant  damage  path  is  clearly  seen  in  the  boxed  regions.  The  damage  size 
progressively  diminishes  with  increasing  sections,  indicating  the  end  of  the  cracked  particles.  The 
particle  area  fraction  (AF),  total  number  of  particles  (NP)  and  total  number  of  cracked  particles 
(NCP)  for  each  section  micrograph  are  presented  in  table  2.  Generally  speaiking,  sections  with  large 
AF  and  NCP  are  found  to  contain  the  larger  cracks.  The  3D  image  by  assembling  2D  micrographs 
in  figure  5.3a  also  shows  the  dominant  damage  path  in  the  boxed  region. 

From  the  microscopic  observation  results,  it  is  found  that  for  the  naturally  aged  material,  the 
main  mode  of  damage  is  by  particle  cracking.  Large  particles  in  particle  rich  regions  are  more 
susceptible  to  cracking  than  those  in  particle  sparse  regions.  Microcracks  in  the  particle  rich  areas 
link  up  to  form  paths  of  dominant  damage.  The  linkage  and  evolution  of  these  larger  cracks  lead 
to  the  overall  failure  of  the  material.  These  paths  are  approximately  perpendicular  to  the  tensile 
loading  direction.  Thus,  spatial  distribution  of  particles  plays  a  more  important  role  in  damage 
than  particle  size  for  this  material. 
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5.3  Equivalent  Microstructure  &:  Mesh  Generation 

The  actual  3D  geometry  of  particles,  as  seen  in  figure  5.3a,  can  be  quite  complex  and  an  exhaustive 
database  is  required  to  store  all  geometric  details.  To  avert  this,  equivalent  microstructures  that 
closely  approximate  the  actual  morphology  but  are  computationally  less  demanding,  are  generated. 
In  this  process,  each  particle  and  an  associated  crack  are  replaced  by  equivalent  ellipses  (in  2D)  or 
ellipsoids  (in  3D).  This  method  economizes  the  image  analysis  and  characterization  process  by  way 
of  well  known  geometric  properties.  For  obtaining  equivalent  microstructures,  digitized  image  data 
is  first  transferred  into  a  binary  format  to  distinguish  between  the  particle,  matrix  and  crack  phases. 
The  0th  (/o),  Tst  {IxJy)  and  2nd  [IxxJyy]  order  geometric  moments  are  then  computed  for  each 
particle  by  adding  contributions  from  each  voxel  (in  3D)  or  pixel  (in  2D)  that  lies  within  the  particle 
boundary.  For  2D  microstructures  the  computed  moments  are  equated  to  the  moment  formulae 
for  ellipses  to  evaluate  the  centroidal  coordinates  (xc,yc)i  half  major  and  minor  axis  lengths  (a,  6), 
and  orientation  0  of  the  major  ajcis  from: 

Xc  =  ~  ~  y 2  'h.sjc'l  -  iC^)  ,  b  =  ^-(C'l  -  -  4C^) 

"  =  <=■') 

where  Ci  =  4  *  and  C2  =  ^.  For  3D  microstructures,  the  centroidal  coordinates 

(xcj/c.Zc)  of  the  equivalent  ellipsoid  are  first  evaluated  from  the  0th  and  1st  order  moments  as: 

=z  ^  ,  j/c  =  ^  principal  directions  (or  orientations  of  the  three  axes)  for  the 

ellipsoids  are  obtained  from  the  eigen-values  of  the  2nd  order  moments  lij  {i  =  1..3,;  =  1..3).  The 
major  (2a),  intermediate  (26)  and  minor  (2c)  axes  of  the  equivalent  ellipsoids  are  then  obtained 
from  the  principal  moments  Ii,l2^h 

a  =  ^ ^(/2  +  I3-  h)  ,  6  =  ^ ^(/i  +  h  -  h)  1  c  =  ^  — (/i  -h  /2  -  h)  (5-2) 

A  simulated  3D  microstructure  with  particles  (grey)  and  cracks  (black)  is  shown  in  figure  5.3b. 
The  microstructures  are  then  tessellated  into  a  mesh  of  2D  and  3D  Voronoi  cells,  by  surface  based 
algorithms  detailed  in  [25,  26].  In  figure  5.3c,  the  mesh  of  Voronoi  cells  is  created  based  on  the 
morphology  of  particles,  while  in  figure  5.3d  the  mesh  is  due  to  tessellation  based  on  the  geometry 
of  particle  cracks.  Tessellation  into  a  mesh  of  Voronoi  coll.s  plays  an  important  role  in  developing 
geometric  descriptors  for  quantitative  characterization.  They  represent  regions  of  immediate  influ¬ 
ence  for  each  heterogeneity  and  also  define  neighbor  of  each  lieterogeneity  from  individual  faces  or 
edges  of  the  Voronoi  cells.  This  facilitates  easy  evaluation  of  parameters  like  local  area  fractions, 
near  neighbor  and  nearest  neighbor  distances  and  orientations. 

5.4  Microstructure  and  Damage  Characterization  '=• 

The  morphology  of  pau-ticles  and  associated  damage  or  microcracks  can  be  characterized  by  various 
functions  of  size,  shape,  orientation  and  spatial  distributi(jn.  .\  number  of  these  classifier  functions 
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have  been  used  by  the  authors  and  others  in  [24,  25,  26,  34,  35,  20,  21,  22,  6]  to  characterize  various 
aspects  of  microstructural  morphology.  In  this  section,  some  of  these  functions  are  considered 
for  the  3D  microstructure  and  2D  micrographs  to  investigate  the  relation  between  morphological 
characteristics  and  the  path  of  dominant  damage  in  the  material.  The  specimen  t6  with  a  large 
microcrack  is  considered  for  this  study. 

In  the  first  exercise,  a  sensitivity  analysis  is  done  with  the  simulated  3D  microstructure  in  figure 
5.3b  to  reveal  the  dependence  of  damage  on  microstructural  variables.  Two  damage  parameters, 
viz.  the  number  fraction  of  cracked  particles  (nf)  and  the  volume  fraction  of  cracked  particles  (vf) 
are  chosen  to  manifest  the  damage  level  in  the  DR.A.  Six  microstructural  parameters  are  considered, 
viz.  (i)  particle  equivalent  size  (diameter);  (ii)  nearest  neighbor  distance  computed  as  the  distance 
between  particles  that  share  a  common  Voronoi  cell  edge;  (iii)  local  volume  fraction  measured  as 
a  ratio  of  the  particle  size  to  that  of  the  associated  Voronoi  cell;  (iv)  particle  shape  or  ellipsoid 
aspect  ratio;  (v)  nearest  neighbor  orientation,  measured  as  the  angle  between  a  line  joining  the 
centers  of  particle  and  its  nearest  neighbor,  and  the  loading  direction;  and  (vi)  particle  orientation 
with  respect  to  the  loading  direction.  The  cracked  particle  fractions  are  plotted  as  functions  of 
these  parameters  in  figure  5.4.  A  linear  interpolation,  obtained  by  a  least  square  fit,  yields  the 
corresponding  overall  gradient  or  slope. 

While  both  the  nf  and  vf  plots  coincide  for  the  particle  size  plot  (i),  large  differences  are  noted 
for  nearest  neighbor  distances  (ii)  and  aspect  ratios  (iv).  Largest  slopes  of  these  plots  are  observed 
with  particle  size,  nearest  neighbor  distance  and  local  volume  fraction.  This  infers  that  the  strongest 
influence  on  particle  cracking  comes  from  the  size  and  local  spatial  distribution.  Particle  shape  has 
a  relatively  smaller  effect  on  damage  initiation.  Sensitivity  of  damage  to  particle  orientation  and 
nearest  neighbor  orientation  is  found  to  be  minimal  for  this  material. 

The  characteristics  of  particles  forming  the  dominant  damage  path  (within  the  marked  box  in 
5.3a)  are  compared  with  those  for  all  cracked  particles  in  the  histograms  of  figure  5.5.  The  dotted 
lines  correspond  to  all  cracked  particles  while  the  shaded  areas  are  for  cracked  particles  in  the 
dominant  damage  region  only.  The  histograms  are  with  respect  to  three  variables  that  are  found 
to  play  important  roles  in  the  damage  process,  viz.  the  particle  size,  nearest  neighbor  distance 
and  orientation  with  respect  to  the  loading  direction.  While  the  range  of  sizes  for  all  cracked 
particles  is  4  ~  I3nm,  that  for  the  particles  forming  the  dominant  damage  path  is  5.7  ~ 

This  reveals  that  larger  particles  generally  contribute  to  dominant  damage  path.  The  plot  for 
nearest  neighbor  clearly  exhibits  the  influence  of  particle  rich  areas  (clustering  or  alignment)  on 
the  preferential  growth  of  damage.  The  nearest  neighbor  distance  for  particles  in  the  dominant 
damage  path  are  in  the  range  0.4  ~  3.7/im  when  compared  with  the  range  0.4  ~  12.3/im  for  all 
cracked  particles.  The  histogr2un  of  cracked  particles  as  a  function  of  the  orientation  with  respect  to 
the  loading  direction  reveals  that  particles  with  major  axis  along  the  loading  direction  (0°  and  180°) 
are  generally  susceptible  to  cracking.  This  is  much  more  pervasive  for  particles  in  the  dominant 
damage  path,  due  to  the  smaller  cross-sectional  areas  normal  to  loading.  In  conclusion,  particles 
in  the  dominamt  damage  path  generally  have  larger  size,  are  in  particle  rich  areas,  and  are  oriented 
in  loading  direction. 

Finally,  it  is  of  interest  to  identify  discriminating  characteristics  of  2D  micrographs  that  may 
be  helpful  in  making  dominant  damage  predictions  for  the  actual  3D  microstructures.  Two  repre- 
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sentative  micrographs,  viz.  section  1  which  contains  a  dominant  damage  and  section  14  with¬ 
out  any  dominant  damage,  but  only  scattered  particle  cracks  are  compared  with  the  3D  mi¬ 
croregion.  Four  important  characterization  functions  viz.  (a)  the  probability  density  function 
of  particle  equivalent  size  (diameter),  (b)  the  probability  density  function  of  the  nearest  neighbor 
distance,  (c)  the  probability  density  function  of  the  local  area/ volume  fraction  and  (d)  a  trans¬ 
formation  function  L(r)  of  a  second  order  intensity  function  K{r),  are  plotted  in  figure  5.6  for 
2D  and  3D  micrographs.  The  second  order  intensity  function  K{r)  and  its  transformed  functions 
in  2D,  and  L(r)  =  in  3D)  capture  second  order  statistics  of  spatial 

distributions  are  used  as  a  graphical  tool  for  detecting  departures  from  a  homogeneous  Poisson 
process  [34.  35,  25,  26].  The  plot  of  L{r)  vs.  r  is  a  45°  straight  line  for  a  pure  Poisson  distribution. 

The  plots  distinctly  reveal  a  few  important  features  of  the  micrographs.  The  particle  size 
distribution  for  the  two  2D  micrographs  are  similar  and  the  tails  are  significantly  shorter  than 
3D.  As  is  expected,  3D  particle  sizes  are  larger  than  2D  particle  section  sizes  due  to  sectioning 
along  non-principal  planes.  However  the  probabilities  of  both  the  nearest  neighbor  distances  and 
local  area  fractions  in  figures  5.6b  and  c  yield  a  distinguishing  characteristic.  The  micrograph  with 
dominant  damage  has  peaks  and  valleys,  as  well  as  tails  that  are  very  similar  to  that  for  3D.  The 
peaks  which  reflect  particle  rich  regions  and  the  tails  which  reflect  sparse  areas  are  both  found  to 
be  important  discriminants.  Deviation  from  the  L{r )  =  r  function  or  the  45°  line  represents  a 
bias  towards  clustering.  The  section  with  the  dominant  damage  has  a  larger  deviation  from  the 
random  distribution  in  comparison  with  the  section  without  major  cracking,  and  is  closer  to  the 
3D  response.  In  summary,  it  may  be  concluded  that  when  analyzing  2D  sections,  the  liklehood  of 
better  representation  of  dominant  damage  are  for  those  sections  that  have  higher  peaks  at  lower 
near  neighbor  distances  with  longer  tails  and  have  higher  deviation  from  the  Poisson  distribution. 
Similar  observations  have  also  been  made  in  [5,  1,  38,  29]. 


5.5  Damage  Simulation  by  Voronoi  Cell  FEM 

Two  dimensional  plane  strain/stress  simulations  of  the  microstructural  damage  evolution  is  con¬ 
ducted  by  the  Voronoi  cell  finite  element  model  (VCFEM)  described  in  [9,  10,  34,  35].  The  current 
2D  VCFEM  only  accommodates  particle  cracking,  and  hence  matrix  cracking  is  ignored  in  the  sim¬ 
ulations.  The  simulations  are  useful  in  understanding  the  damage  evolution  process  by  a  sequence 
of  particle  cracking.  Rectangular  195/im  x  155.018/im  micrographs  as  shown  in  figure  5.8a,b  are 
analyzed  with  monotonically  increasing  strains.  Periodicity  boundary  conditions  are  imposed  by 
requiring  edges  to  remain  straight  and  parallel  to  the  original  direction  throughout  deformation  as: 

Ux  =  0  {on  1  =  0),  Uy  =  0  {on  y  =  0)  ,  Ui  =  n^p  {on  x  =  Lx)  ,  Uy  =  Dy  {on  y  =  Ly) 

Ty  =  0  {on  X  =  0/Lx)  ,  Tx  =  0  {on  y  =  0/Ly)  (5-3) 

where  Uap  is  an  applied  displacement  and  D*  is  determined  from  the  average  force  condition 
/ .  Txdx  =  0  on  y  =  Ly.  The  reinforcing  phase  of  SiC  particles  are  assumed  ts^be  brittle  and 
is  modeled  with  the  linear  elastic  properties:  Young  s  modulus  E  =  427  GPa,  Poisson’s  ratio  i/ 
=  0.17.  The  aluminum  matrix  material  is  assumed  to  be  ductile  and  is  modeled  by  small  defor¬ 
mation  isotropic  hardening  J2  elasto-plasticity  theory  with  properties:  Young’s  modulus  E  =  72 
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GPa,  Poisson’s  ratio  u  =  0.33  and  the  post  yield  elastic-plastic  behavior  is  obtained  from  [39]  as 
shown  in  Bgure  5.7.  Microstructural  damage  by  particle  cracking  is  assumed  to  be  governed  by 
a  maximum  principal  stress  or  Rankine  criterion.  In  this  criterion,  a  crack  is  initiated  when  the 
maximum  principal  stress  in  tension  exceeds  a  critical  fracture  stress  Go-  at  a  point.  The  crack  is 
oriented  at  right  angle  to  the  principal  stress  direction.  The  critical  stress  Gct  is  also  influenced  by 
the  particle  size  due  to  the  existence  of  microcracks.  To  account  for  the  size  effect  in  ctct,  a  Weibull 
distribution  based  criterion  is  used,  where  the  probability  of  particle  fracture  Pf{A,  a)  is  related  to 
the  particle  volume/area  v  and  the  maximum  principal  stress  a/  as: 

Pf{v,aj)  =  1  -  eip[-u(  — )'”]  (5.4) 

where  gq  and  m  are  two  material  parameters  in  the  Weibull  distribution  that  are  calibrated  from 
experiments. 

5.5.1  Calibration  of  Weibull  parameters  gq  and  m 

In  the  two  parameter  Weibull  model,  the  fraction  of  fractured  particles  may  be  obtained  (see 
[41.  42.  26])  from  a  known  probability  distribution  of  particle  volumes  p{v),  as: 

p{v)=  l^p{v)Pf(v,Gi)dv  ss  j^p(t;,)  (1  -  exp[-^(^)'"])  Avj  (5.5) 

where  p(vi)  is  the  probability  density  distribution  of  particle  volume/area  Uj.  The  entire  area  is 
divided  into  N  intervals  such  that  Auj  =  u,  -  Uj-i,  g\  is  the  average  particle  maximum  principal 
stress  for  particles  with  size  in  the  range  of  [uj-i,  Uj]  and  uq  is  a  reference  area  taken  to  be  the  average 
area.  The  fraction  of  cracked  particles  p  is  readily  obtained  from  the  experimental  micrographs. 
.4gain,  the  section  micrographs  2,  8  and  14  are  used  to  calibrate  the  Weibull  parameters.  The 
fractions  of  cracked  particles  and  the  average  pau-ticle  area  for  these  three  sections  cire  31.78%, 
24.76%,  28.57%  and  53.43,  48.91  and  52.67  pvn?  respectively.  The  maximum  principal  stress  g)  for 
each  particle  is  obtained  from  VCFEM  simulation  prior  to  the  onset  of  particle  cracking  at  a  true 
strain  of  c  =  8.88%.  Prom  the  experimental  observations  it  is  assumed  that  no  major  damage  has 
initiated  at  this  strain.  The  Weibull  parameter  m  is  assumed  to  take  integer  values  between  1  and 
8  following  [42.  26]  and  the  corresponding  values  of  gq  are  given  in  table  3. 

The  Weibull  pjirameters  are  also  calibrated  using  a  3D  .ABAQllS  model  simulation  of  a  cubic 
unit  cell  with  a  single,  15%  volume  fraction,  spherical  particle  as  described  in  [26].  The  1x1x1  unit 
cell  model  has  a  particle  of  radius  R  =  0.66.  A  modified  form  of  equation  5.5  is  used  to  account 
for  the  shape  veiriability  of  the  particles  as 

PCimar  r^rnax 

p{a,v)=  /  p(a)p{v)  Pf{v.Gi)dv  (5.6) 

(^m%n  ^min  ^ 

where  a  corresponds  to  the  particle  aspect  ratio.  The  particle  size  and  shape  distribution  functions 
p(v)  and  p(a)  axe  calculated  from  the  computer  simulated  representation  of  the  actual  3D  mi¬ 
crostructure  shown  in  figure  5.3.  This  average  particle  volume  and  the  fraction  of  cracked  particle 
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are  directly  computed  from  figure  5.3  as  v  =642. and  p  =  45.48%.  The  average  particle  stress 
at  a  macroscopic  strain  e  =  8.88%  is  obtained  from  the  ABAQUS  simulation  as  cXp  =  862.60  MPa. 
Results  of  calibration  with  and  without  shape  effects  cure  documented  in  table  3.  It  is  found  that  the 
best  agreement  in  ctq  for  all  2D  sections  and  3D  is  obtained  for  m  between  4  and  5.  Consequently 
the  parameter  is  chosen  to  be  m  =  4.2.  The  corresponding  value  of  oq  for  section  2  is  3.04  GPa,  for 
section  8  is  3.19  GPa,  for  section  14  is  2.79  GPa  and  the  average  of  these  sections  is  ao  =  Z.OlGPa. 

Results  of  VCFEM  analysis  of  the  simulated  micrographs  of  sections  1.  3,  5  and  9  are  provided 
in  table  4.  The  number  of  cracked  particles  at  a  macroscopic  strain  of  8.88%  by  VCFEM  are 
compared  with  experimental  results.  While  the  general  agreement  is  quite  good,  it  is  seen  that  the 
concurrence  is  particularly  favorable  when  the  simulation  is  conducted  with  a  ctq  that  is  obtained 
from  a  section  that  is  near  to  the  one  being  analyzed.  For  example,  the  results  of  sections  1  and 
3  are  very  good  when  ao  =  3.01  GPa,  which  is  obtained  from  section  2.  This  concurrence  may 
be  attributed  to  the  similarity  in  the  distribution  of  heterogeneities  in  neighboring  sections,  and 
suggests  that  spatial  distribution  has  a  strong  effect  on  the  Weibull  parameters. 

Microscopic  Damage  Analysis 

Various  results  for  section  1  which  contain  a  dominant  damage  path  are  generated  by  VCFEM 
simulation  and  compared  with  experimental  observations  in  this  section.  The  macroscopic  stress- 
strain  plot  for  plane  strain  and  plane  stress  assumptions  are  compared  with  experimental  results  in 
figure  5.7a.  The  overall  yield  strength  is  better  predicted  by  the  plane  stress  model.  However,  the 
post  yield  behavior  with  plane  strain  conditions  is  much  closer  to  the  experimental  results.  The 
initial  higher  yield  strength  is  expected  with  plane  strain  due  to  the  plastic  constraint  caused  by 
the  =  0  condition.  A  shifted  stress-strain  plot  (modified  plane  strain  VCFEM  result  in  figure 
5.7a)  where  the  stresses  are  reduced  by  the  initial  difference  in  yield  stress  shows  a  very  good 
match  between  experiments  and  simulation.  Thus  plane  strain  assumptions  are  used  in  subsequent 
computations.  Figure  5.7b  is  intended  to  predict  the  onset  of  plastic  instability  by  the  model  and 
compare  it  with  actual  fracture  observed  in  the  experiments.  The  use  of  the  Considere  criterion 
to  predict  the  onset  of  plastic  instability  has  been  suggested  by  Llorca  [2,  41]  in  the  absence  of 
dilatational  strain  associated  with  reinforcement  fracture.  In  this  criterion,  the  average  stress  a  is 
related  to  the  strain  hardening  rate  ^  as 

a  =  ^  (5.7) 

de 

The  strain  derived  from  this  relation  corresponds  to  the  lower  bound  of  the  tensile  ductility  since 
it  controls  the  composite  load  bearing  capacity.  Three  sets  of  curves  are  plotted  in  the  figure  5.7b 
corresponding  to  the  matrix  material,  the  VCFEM  results  in  plane  strain  and  the  experimental 
results.  It  is  seen  that  the  Considere  criterion  (junction  of  the  two  curves)  predicts  the  experimental 
point  corresponding  to  the  onset  of  fracture  rather  well.  Additionally  the  2D  prediction  of  the  plane 
strain  simulation  is  adso  quite  good  and  can  be  used  with  reasonable  confidence.  ^ 

Microstructural  results  of  the  simulation  are  compared  with  experiments  in  figure  5.8.  The 
computed  micrograph  with  evolved  damage  for  section  1  is  compared  with  the  experimental  micro¬ 
graph  at  8.88%  strain  in  figure  5.8a  and  b.  The  damaged  particles  are  shown  with  the  contained 
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crack.  Most  damaged  particles  in  the  simulation  coincide  with  the  experimental  results,  with  the 
box  indicating  the  dominant  damage  path.  The  damage  path  is  approximately  perpendicular  to 
the  tensile  loading  direction.  Figure  5.8b  also  shows  the  contour  plot  of  eflFective  plastic  strain 
in  the  ductile  matrix  and  indicates  the  path  of  damage  linkage.  The  plastic  strain  is  higher  and 
localized  between  cracked  particles  and  this  is  expected  to  cause  matrbc  cracking.  The  number 
fraction  of  cracked  particles  in  different  size  ranges  are  plotted  in  figure  5.8c.  Very  good  agreement 
is  seen  between  simulation  and  experimental  results.  Figure  5.8d  is  a  contour  plot  of  the  particle 
fracture  probability  at  8.88%  strain.  The  black  shade  corresponds  to  the  highest  probability  and 
fractured  paiticles  are  illustrated  in  white  with  a  crack.  Similar  plots  (not  shown)  at  earlier  stages 
of  deformation  show  that  several  particles  with  higher  probability  at  the  smaller  strain  have  cracked 
with  deformation.  The  number  fraction  of  cracked  particles  as  a  function  of  straining  are  plotted 
for  sections  1,  5  and  9  together  with  the  experimental  observation  in  figure  5.9.  At  lower  strains  the 
number  fraction  of  of  cracked  particles  for  sections  1  and  5  with  particle  rich  regions  are  higher  than 
that  for  section  9.  This  is  due  to  higher  stress  concentrations  particle  rich  areas  that  are  enough  to 
fracture  some  particles  even  at  low  strains.  With  increasing  strains  more  particles  start  to  crack  in 
the  section  9  and  exceeds  that  for  section  5  which  has  less  particles  in  the  clustered  regions.  The 
2D  simulations  however  exhibit  less  cracked  particles  than  that  in  the  actual  3D  microstructure. 


5.6  Characteristic  Size  of  Microstructures 

The  influence  region  of  local  morphology  on  the  mechanical  response  is  characterized  by  a  mi- 
crostructural  representative  material  element  (RME)  that  is  critical  in  delineating  length  scales. 
The  RME  depicts  a  region  which  is  assumed  to  be  representative  of  the  entire  microstructufe. 
Functions  that  distinguish  between  variations  in  stress/strain  distributions  for  local  disturbances 
in  microstructural  patterns  can  provide  important  insight  on  microstructure-property  relations. 
Marked  correlation  functions,  discussed  in  [22,  10,  35,  26]  for  multivariate  characterization  ofpat- 
terns,  are  evaluated  to  characterize  length  scales  or  RME  size  in  the  presence  of  damage.  A  mark 
may  be  identified  with  an  appropriate  microstructural  variables,  e.g.  in  this  case  a  variable  that 
related  to  quantification  of  damage.  The  marked  correlation  function  for  a  heterogeneous  domain 
W  of  volume  V  containing  N  heterogeneities  is  mathematically  expressed  as  [22,  23,  10): 


M(r)  = 


l^l/(4«r^) 


/f(r)  =  and  g{r)  = 


1=1 k=l 


1  dKjr) 
4irr'^  dr 


(5.8) 


where  K{r)  is  the  second  order  intensity  function  defined  in  i34,  35],  g{r)  is  the  pair  distribution 
function  and  H(r)  is  the  mark  intensity  function.  The  H{r)  function  reduces  to  the  K{r)  function  if 
all  heterogeneities  have  the  same  mark.  A  mark  associated  with  the  heterogeneity  is  denoted  as 
m,.  k'  is  the  number  of  heterogeneities  which  have  their  (  (Miters  within  a  sphere  of  radius  r  around 
heterogeneity,  for  which  the  mark  is  m/t,  and  m  is  the  mean  of  all  marks.  By  definition  M{r) 
establishes  a  relation  between  the  location  and  associated  variables  for  heterogeneities.  Two  marks 
are  considered  in  this  study.  The  first  corresponds  to  particle  cracks  and  are  designated  as  mi  =  1 
for  a  cracked  particle  and  mj  =  2  for  an  intact  particle.  The  second  corresponds  to  the  probability 
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of  particle  fracture,  which  signifies  the  propensity  to  advance  the  microstructural  damage  state. 
The  M(r)  function  statistically  stabilizes  at  near-unit  values  at  a  distance  Tinter  at  which  the  local 
morphology  ceases  to  have  any  significant  infiuence  on  evolving  variable.  Values  of  M(r)  >  1  show 
positive  correlation,  while  M{r)  <  1  indicates  repulsion  between  marks.  This  distance  rinter  is  an 
indication  of  the  physical  range  of  interaction  and  is  significant  in  making  decisions  about  length 
scales  and  the  RME  size. 

The  marked  correlation  functions  corresponding  to  cracked  particles  and  probability  of  fracture, 
are  plotted  in  figures  5.11a  and  b  from  the  simulation  of  the  entire  micrograph  of  section  1  (termed 
as  RME  0)  with  dimensions  195/'7r<.  x  The  dotted  line  corresponds  to  the  unit  M(r) 

for  uniform  distribution  of  spherical  heterogeneities  with  identical  marks.  Contour  plots  of  the 
equivalent  plastic  strain  in  the  simulated  micrograph  with  cracked  particles  are  shown  in  figure 
5.10a.  The  particle  area  fraction  for  this  micrograph  18.37%  and  the  total  number  of  particles 
and  cracked  particles  are  105  and  34  respectively.  The  plots  are  made  with  only  upto  40%  of 
the  entire  micrograph,  or  SO/xm  to  avoid  boundary  effects  in  M(r).  The  M(r)  functions  in  both 
figures  approximately  stabilize  at  near-unit  values  at  a  distance  Tinter  of  about  60  At  this 
distance,  the  local  morphology  is  expected  to  have  a  significantly  reduced  influence  on  the  evolving 
variables.  The  slower  attenuation  of  M(r)  for  particle  fracture  at  shorter  range  indicates  the  strong 
effect  of  the  local  morphology  on  damage  evolution.  Next,  a  smaller  region  (RME  1)  is  selected 
for  damage  simulation  corresponding  to  the  stabilized  region  in  the  M(r)  plots.  Since  the  stable 
region  is  60/im,  the  dimension  of  the  micrograph  is  chosen  to  be  ISO/xm  x  155^m.  incorporating 
the  scaling  factor,  i.e.  ^  =  150.  This  is  shown  with  the  box  in  figure  5.10a.  Again  the  contour 
plots  of  plastic  strain  with  cracked  particles  by  VCFEM  simulation  are  shown  in  figure  5.10b..  The 
dominant  crack  behavior  is  quite  similar  to  that  for  RME  0,  even  though  there  is  some  difference 
near  the  boundary.  Also  the  plastic  strain  contours  and  limiting  values  are  similar.  The  particle 
area  fraction  for  RME  1  is  18.13%  with  the  total  number  of  particles  and  cracked  particles  at  84 
and  28  respectively.  The  Af(r)  plots  in  figure  5.11c  and  d  show  that  the  functions  may  still  be 
assumed  to  stabilize  at  around  fiO^m.  A  smaller  subset  (RME  2),  with  dimensions  116MTnxll5/im 
is  next  simulated  and  the  plastic  strain  is  depicted  in  figure  5.10c.  Significantly  different  plastic 
strains  and  cracking  pattern  is  observed  for  this  microstructure.  The  particle  area  fraction  for 
this  micrograph  is  18.68%  with  a  total  of  51  particles  of  which  18  are  cracked.  The  plots  of  M(r) 
function  do  not  stabilize  in  the  domain  of  the  simulation  window.  Through  this  analysis  the  size 
effect  of  microstructure,  needed  for  adequate  representation  and  analysis  in  the  presence  of  evolving 
damage  is  demonstrated. 

5.7  Conclusions 

In  this  work,  a  combination  of  experimental  and  computational  methods  are  utilized  to  characterize 
and  understand  the  evolution  of  microscopic  damage  that  cause  failure  in  naturally  aged  commercial 
SiC  particle  reinforced  DRA’s.  The  main  mode  of  damage  for  the  naturally  aged  gaterial  is  found 
to  be  particle  cracking.  Larger  particles  in  particle  rich  regions  are  more  susceptible  to  cracking 
than  those  in  particle  spairse  regions.  Spatial  distribution  of  particles  plays  a  more  important  role 
in  damage  than  particle  size  for  this  material.  A  sensitivity  analysis  with  respect  to  microstructural 
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parameters  infers  that  the  strongest  influence  on  particle  cracking  comes  from  the  size  and  local 
spatial  distribution.  Particle  shape,  orientation  and  nearest  neighbor  orientation  have  relatively 
smaller  effect  on  damage  initiation.  Histograms  of  particles  forming  the  dominant  damage  path  in 
comparison  with  all  cracked  particles  reveal  that  larger  particles  oriented  in  loading  direction  and  in 
relatively  rich  areas  are  more  susceptible  to  contribute  to  a  dominant  crack  in  the  microstructure. 
In  an  attempt  to  identify  discriminating  characteristics  of  2D  micrographs  that  may  be  of  helpful 
in  making  dominant  damage  predictions  for  the  actual  3D  microstructures,  probability  density 
functions  of  particle  size  nearest  neighbor  distance,  and  second  order  intensity  function  K{r)  of 
spatial  distribution  are  plotted.  Better  representation  of  damage  is  possible  with  those  sections 
that  have  higher  peaks  at  lower  near  neighbor  distances  and  longer  tails,  as  well  as  have  propensity 
towards  clustering. 

Next  the  two  dimensional  Voronoi  cell  finite  element  model  is  used  to  simulate  microstructural 
damage  evolution  in  computer  generated  equivalent  micrographs.  Both  macroscopic  and  micro¬ 
scopic  variables  obtained  by  the  VCFEM  simulation  are  compared  with  experimental  observations. 
The  macroscopic  stress-strain  plot  for  the  plane  strain  analysis  is  found  to  yield  quite  good  match 
with  experiments  if  the  difference  in  the  initial  yield  strength  due  to  plastic  constraint  is  subtracted. 
Prediction  of  the  onset  of  plastic  instability  by  the  Considere  criterion  is  also  found  to  be  in  reason¬ 
ably  good  agreement  with  the  experimental  results.  For  the  microstructural  results  with  number 
of  cracked  particles  in  different  size  ranges,  the  Weibull  model  is  found  to  give. better  concurrence 
with  experiments.  A  plot  of  the  number  fraction  of  cracked  particles  as  a  function  of  straining 
shows  that  at  lower  strains  sections  with  particle  rich  regions  damage  rapidly,  but  the  rate  slows 
down  with  additional  deformation.  Finally,  the  marked  correlation  functions  are  evaluated  to  char¬ 
acterize  length  scales  and  representative  material  element  size  in  the  presence  of  damage.  Particle 
cracks  and  the  probability  of  particle  fracture  are  chosen  to  be  the  marks.  The  study  reveals  that  a 
significantly  large  portion  of  the  microstructure  should  be  analyzed  for  reasonable  accuracy  in  the 
presence  of  damage.  The  correlation  functions  do  not  stabilize  below  a  certain  length  scales  and 
this  keeps  growing  with  increased  damage.  In  summary,  various  important  characteristics  of  grow¬ 
ing  damage  are  investigated  in  this  work  to  understand  the  role  of  microstructure  in  the  material 
failure  process. 
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Figure  5.1:  (a)  Interrupted  uniaxial  tensile  test  specimen  for  naturally  aged  material  and  sample 
coupon  for  serial  sectioning(unit  in  mm);  (b)  Load-strain  plots  for  two  specimens  of  naturally  aged 
DRA.  Dark  points  indicate  where  the  loading  is  interrupted  or  where  the  specimen  is  fractured. 
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Figure  5.2;  Micrographs  of  different  sections  of  the  t6  specimen  showing  cracked  particles;  (a) 
section  3,  (b)  section  5. 
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Figure  5.3:  3D  microstructure  for  SiC  particle  reinforced  DRA  after  interrupted  test  (units  in 
^m):  (a)  A  computer  created  image  by  seriallj;St3.cking  section  micrographs,  (b)  Simulated  mi¬ 
crostructure  of  ellipsoidal  particles  and  cracks.  Tessellation  for  2D  and  3D  microstructures,  (c)  3D 
tessellation  based  on  particle  morphology  and  (d)  3D  tessellation  based  on  microcrack  morphology. 
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Figure  5.4:  Sensitivity  of  damage  to  various  microstructural  variables 
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Figure  5.5:  Histograms  comparing  characteristics  of  particles  with  dominant  damage  with  all 
cracked  particles 
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Figxire  5.6;  Characterization  functions  for  2D  sections  and  3D  microstructure. 
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Figure  5.7:  (a)  Macroscopic  stress-strain  response  by  plane  strain  and  plane  stress  VCFEM  simu¬ 
lation,  (b)  Stress-strzun  hardening  rate  plots  for  the  Considere  condition. 
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Figure  5.8:  (a)  Experimental  micrographs,  (b)  VCFEM  simulated  micrograph  showing  damage  and 
contour  plot  of  effective  plastic  strain  at  8.88%  strain  in  section  1,  (c)  histogram  oHiumber  fraction 
of  craudeed  pzirticles  as  a  function  of  p2irticle  size  by  Weibull  based  probabilistic  criterion,  and  (d) 
contour  plot  of  particle  fracture  probability  of  section  1  at  8.88%  strain. 
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Figure  5.9:  Number  fraction  of  cracked  particles  as  a  function  of  straining. 
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Figure  5.10:  Simulations  showing  effective  plastic  strain(%)  and  cracked  particles  in  three  different 
subsets  of  the  entire  micrograph  of  section  1,  (a)  RME  0  with  dimension  195/iTn  x‘l-55/xTn  (b)  RME 
1  with  dimension  150/im  x  155/im  and  (c)  RME  2  with  dimension  116/xTnxll5MTn. 
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Figure  5.11:  Marked  correlation  function  as  a  function  of  radial  distance;  (a)  cracked  particles  as 
marks  and  (b)  probability  of  cracking  as  makr  for  RME  0;  (c)  cracked  particles  as  marks  and  (d) 
probability  of  cracking  as  makr  for  RME  1;  (e)  cracked  particles  as  marks  and.  (f)  probability  of 
cracking  as  makr  for  RME  2. 
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